new; cls;
/* MC Bootstrap Lookback Call & Put Options            */
/*     Data from Sheldon M. Ross(1999) Table 10.5 in   */
/*     "An Introduction to Mathematical Finance"       */
/* (It will take some time. Wait for a long moment...) */
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
;
data=P;
n=rows(data);
ST=data[n];
v={-0.02,-0.01,0,0.01,0.02};         /* 5 break points */
K=20;
r=0.03;
t=2.9;
times=200;                                 /* Try 5000 */


print/lz "      n=" n "(its result is n-1)";
print/lz "  times=" times;
print "C=" MClookbackC(ST,data,v,r,t,times);
print "P=" MClookbackP(ST,data,v,r,t,times);
print;
print "Fixed Strike version:";
print "C=" MClookbackCf(ST,data,v,K,r,t,times);
print "P=" MClookbackPf(ST,data,v,K,r,t,times);




/*
**  mclookback.txt - MC Bootstrap Lookback Call & Put Options.
** (C) Copyright 2004 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Lookback option prices by MC bootstrapped simulation
**              from real past data.
**
**  Format:     C=MClookbackC(ST,data,v,r,t,times);
**              P=MClookbackP(ST,data,v,r,t,times);
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**             
**              v        vector,  5 x 1 vector of break points
**
**              r        scalar,  risk-free interest rate
**
**              t        scalar,  maturity (the same as data period)
**
**              times    scalar,  times of bootstrap
**
**
**  Output:     C        scalar,  call option price
**              P        scalar,  put  option price
**
**  Notice:     It takes very long time to run. Try some small 'times' first.
**
*/

proc MClookbackC(ST,data,v,r,t,times);
   local n,OV,i,S,C;
   n=rows(data)-1;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=mcsampler66(ST,data,v);
      OV[i]=S[n]-minc(S[2:n]);
      i=i+1;
   endo;
   C=exp(-r*t)*meanc(maxc((OV~zeros(times,1))'));
   retp(C);
endp;

proc MClookbackP(ST,data,v,r,t,times);
   local n,OV,i,S,P;
   n=rows(data)-1;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=mcsampler66(ST,data,v);
      OV[i]=maxc(S[2:n])-S[n];
      i=i+1;
   endo;
   P=exp(-r*t)*meanc(maxc((OV~zeros(times,1))'));
   retp(P);
endp;

/*
**  mclookback.txt - Fixed Strike version of MC Bootstrap Lookback Call & Put Options.
** (C) Copyright 2004 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Lookback option prices by MC bootstrapped simulation
**              from real past data.
**
**  Format:     C=MClookbackCf(ST,data,v,K,r,t,times);
**              P=MClookbackPf(ST,data,v,K,r,t,times);
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**             
**              v        vector,  5 x 1 vector of break points
**
**              K        scalar,  strike price
**
**              r        scalar,  risk-free interest rate
**
**              t        scalar,  maturity (the same as data period)
**
**              times    scalar,  times of bootstrap
**
**
**  Output:     C        scalar,  call option price
**              P        scalar,  put  option price
**
**  Notice:     It takes very long time to run. Try some small 'times' first.
**
*/

proc MClookbackCf(ST,data,v,K,r,t,times);
   local n,OV,i,S,C;
   n=rows(data)-1;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=mcsampler66(ST,data,v);
      OV[i]=maxc(S[2:n])-K;
      i=i+1;
   endo;
   C=exp(-r*t)*meanc(maxc((OV~zeros(times,1))'));
   retp(C);
endp;

proc MClookbackPf(ST,data,v,K,r,t,times);
   local n,OV,i,S,P;
   n=rows(data)-1;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=mcsampler66(ST,data,v);
      OV[i]=K-minc(S[2:n]);
      i=i+1;
   endo;
   P=exp(-r*t)*meanc(maxc((OV~zeros(times,1))'));
   retp(P);
endp;

/*
**  mcsampler.txt - Markov chain(discrete 4x4 Gibbs sampling) Sampler from real data.
** (C) Copyright 2004 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Samples just one MC bootstrap representative data from real data.
**
**  Format:     S=mcsampler44(ST,data,v);
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**
**              v        vector,  3 x 1 vector of break points
**             
**
**  Output:     S        vector,  n-1 x 1 vector of MC bootstrapped representative data
**
*/

proc mcsampler44(ST,data,v);
   local D,n,DI,i,index,temp,m,joint,a,b,r,S;
   local x11,x12,x13,x14,x21,x22,x23,x24,x31,x32,x33,x34,x41,x42,x43,x44;
   D=ln(data[2:rows(data)])-ln(data[1:rows(data)-1]);
   n=rows(D)-1;
/* Put dummy numbers to each data depending on its state. */
   DI=zeros(n,1);
   i=1;
   do while i<=n;
      if D[i]<=v[1];
         if D[i+1]<=v[1];
            DI[i]=11;
         elseif v[1]<D[i+1] and D[i+1]<=v[2];
            DI[i]=12;
         elseif v[2]<D[i+1] and D[i+1]<=v[3];
            DI[i]=13;
         else;
            DI[i]=14;
         endif;
      elseif v[1]<D[i] and D[i]<=v[2];
         if D[i+1]<=v[1];
            DI[i]=21;
         elseif v[1]<D[i+1] and D[i+1]<=v[2];
            DI[i]=22;
         elseif v[2]<D[i+1] and D[i+1]<=v[3];
            DI[i]=23;
         else;
            DI[i]=24;
         endif;
      elseif v[2]<D[i] and D[i]<=v[3];
         if D[i+1]<=v[1];
            DI[i]=31;
         elseif v[1]<D[i+1] and D[i+1]<=v[2];
            DI[i]=32;
         elseif v[2]<D[i+1] and D[i+1]<=v[3];
            DI[i]=33;
         else;
            DI[i]=34;
         endif;
      else;
         if D[i+1]<=v[1];
            DI[i]=41;
         elseif v[1]<D[i+1] and D[i+1]<=v[2];
            DI[i]=42;
         elseif v[2]<D[i+1] and D[i+1]<=v[3];
            DI[i]=43;
         else;
            DI[i]=44;
         endif;
      endif;
      i=i+1;
   endo;
/* Make subsets of data for 4x4 states. */
   DI=0|DI;
   x11=selif(D,(DI.==11));
   x12=selif(D,(DI.==12));
   x13=selif(D,(DI.==13));
   x14=selif(D,(DI.==14));
   x21=selif(D,(DI.==21));
   x22=selif(D,(DI.==22));
   x23=selif(D,(DI.==23));
   x24=selif(D,(DI.==24));
   x31=selif(D,(DI.==31));
   x32=selif(D,(DI.==32));
   x33=selif(D,(DI.==33));
   x34=selif(D,(DI.==34));
   x41=selif(D,(DI.==41));
   x42=selif(D,(DI.==42));
   x43=selif(D,(DI.==43));
   x44=selif(D,(DI.==44));
/* joint distribution */
   m=rows(x11)~rows(x12)~rows(x13)~rows(x14)|
     rows(x21)~rows(x22)~rows(x23)~rows(x24)|
     rows(x31)~rows(x32)~rows(x33)~rows(x34)|
     rows(x41)~rows(x42)~rows(x43)~rows(x44);
   joint=m/n;
/* Gibbs sampling */
      a=1; b=1;           /* beginning from (1,1) you could change depending on data. */
      temp=zeros(n,1);
      i=1;
      do while i<=n;
         a=usample(seqa(1,1,rows(joint)),joint[.,b]/sumc(joint[.,b]));
         b=usample(seqa(1,1,cols(joint)),joint[a,.]'/sumc(joint[a,.]'));
         if a==1 and b==1;
            index=ceil(rows(x11)*rndu(1,1));
            temp[i]=x11[index];
         elseif a==1 and b==2;
            index=ceil(rows(x12)*rndu(1,1));
            temp[i]=x12[index];     
         elseif a==1 and b==3;
            index=ceil(rows(x13)*rndu(1,1));
            temp[i]=x13[index];             
         elseif a==1 and b==4;
            index=ceil(rows(x14)*rndu(1,1));
            temp[i]=x14[index];             
         elseif a==2 and b==1;
            index=ceil(rows(x21)*rndu(1,1));
            temp[i]=x21[index];         
         elseif a==2 and b==2;
            index=ceil(rows(x22)*rndu(1,1));
            temp[i]=x22[index];         
         elseif a==2 and b==3;
            index=ceil(rows(x23)*rndu(1,1));
            temp[i]=x23[index];         
         elseif a==2 and b==4;
            index=ceil(rows(x24)*rndu(1,1));
            temp[i]=x24[index];         
         elseif a==3 and b==1;
            index=ceil(rows(x31)*rndu(1,1));
            temp[i]=x31[index];         
         elseif a==3 and b==2;
            index=ceil(rows(x32)*rndu(1,1));
            temp[i]=x32[index];         
         elseif a==3 and b==3;
            index=ceil(rows(x33)*rndu(1,1));
            temp[i]=x33[index];         
         elseif a==3 and b==4;
            index=ceil(rows(x34)*rndu(1,1));
            temp[i]=x34[index];         
         elseif a==4 and b==1;
            index=ceil(rows(x41)*rndu(1,1));
            temp[i]=x41[index]; 
         elseif a==4 and b==2;
            index=ceil(rows(x42)*rndu(1,1));
            temp[i]=x42[index];          
         elseif a==4 and b==3;
            index=ceil(rows(x43)*rndu(1,1));
            temp[i]=x43[index];          
         else;
            index=ceil(rows(x44)*rndu(1,1));
            temp[i]=x44[index]; 
         endif;           
         i=i+1;
      endo;
   r=ones(n,1)+temp;
   S=ST*cumprodc(r);
   retp(ST|S);
endp;

/*
**  mcsampler.txt - Markov chain(discrete 6x6 Gibbs sampling) Sampler from real data.
** (C) Copyright 2004 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Samples just one MC bootstrap representative data from real data.
**
**  Format:     S=mcsampler66(ST,data,v);
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**
**              v        vector,  5 x 1 vector of break points
**             
**
**  Output:     S        vector,  n-1 x 1 vector of MC bootstrapped representative data
**
*/

proc mcsampler66(ST,data,v);
   local D,n,DI,i,index,temp,m,joint,a,b,r,S;
   local x11,x12,x13,x14,x15,x16,x21,x22,x23,x24,x25,x26,x31,x32,x33,x34;
   local x35,x36,x41,x42,x43,x44,x45,x46,x51,x52,x53,x54,x55,x56,x61,x62,x63,x64,x65,x66;
   D=ln(data[2:rows(data)])-ln(data[1:rows(data)-1]);
   n=rows(D)-1;
/* Put dummy numbers to each data depending on its state. */
   DI=zeros(n,1);
   i=1;
   do while i<=n;
      if D[i]<=v[1];
         if D[i+1]<=v[1];
            DI[i]=11;
         elseif v[1]<D[i+1] and D[i+1]<=v[2];
            DI[i]=12;
         elseif v[2]<D[i+1] and D[i+1]<=v[3];
            DI[i]=13;
         elseif v[3]<D[i+1] and D[i+1]<=v[4];
            DI[i]=14;
         elseif v[4]<D[i+1] and D[i+1]<=v[5];
            DI[i]=15;
         else;
            DI[i]=16;
         endif;
      elseif v[1]<D[i] and D[i]<=v[2];
         if D[i+1]<=v[1];
            DI[i]=21;
         elseif v[1]<D[i+1] and D[i+1]<=v[2];
            DI[i]=22;
         elseif v[2]<D[i+1] and D[i+1]<=v[3];
            DI[i]=23;
         elseif v[3]<D[i+1] and D[i+1]<=v[4];
            DI[i]=24;
         elseif v[4]<D[i+1] and D[i+1]<=v[5];
            DI[i]=25;
         else;
            DI[i]=26;
         endif;
      elseif v[2]<D[i] and D[i]<=v[3];
         if D[i+1]<=v[1];
            DI[i]=31;
         elseif v[1]<D[i+1] and D[i+1]<=v[2];
            DI[i]=32;
         elseif v[2]<D[i+1] and D[i+1]<=v[3];
            DI[i]=33;
         elseif v[3]<D[i+1] and D[i+1]<=v[4];
            DI[i]=34;
         elseif v[4]<D[i+1] and D[i+1]<=v[5];
            DI[i]=35;
         else;
            DI[i]=36;
         endif;
      elseif v[3]<D[i] and D[i]<=v[4];
         if D[i+1]<=v[1];
            DI[i]=41;
         elseif v[1]<D[i+1] and D[i+1]<=v[2];
            DI[i]=42;
         elseif v[2]<D[i+1] and D[i+1]<=v[3];
            DI[i]=43;
         elseif v[3]<D[i+1] and D[i+1]<=v[4];
            DI[i]=44;
         elseif v[4]<D[i+1] and D[i+1]<=v[5];
            DI[i]=45;
         else;
            DI[i]=46;
         endif;
      elseif v[4]<D[i] and D[i]<=v[5];
         if D[i+1]<=v[1];
            DI[i]=51;
         elseif v[1]<D[i+1] and D[i+1]<=v[2];
            DI[i]=52;
         elseif v[2]<D[i+1] and D[i+1]<=v[3];
            DI[i]=53;
         elseif v[3]<D[i+1] and D[i+1]<=v[4];
            DI[i]=54;
         elseif v[4]<D[i+1] and D[i+1]<=v[5];
            DI[i]=55;
         else;
            DI[i]=56;
         endif;      
      else;
         if D[i+1]<=v[1];
            DI[i]=61;
         elseif v[1]<D[i+1] and D[i+1]<=v[2];
            DI[i]=62;
         elseif v[2]<D[i+1] and D[i+1]<=v[3];
            DI[i]=63;
         elseif v[3]<D[i+1] and D[i+1]<=v[4];
            DI[i]=64;
         elseif v[4]<D[i+1] and D[i+1]<=v[5];
            DI[i]=65;
         else;
            DI[i]=66;
         endif;
      endif;
      i=i+1;
   endo;
/* Make subsets of data for 6x6 states. */
   DI=0|DI;
   x11=selif(D,(DI.==11));
   x12=selif(D,(DI.==12));
   x13=selif(D,(DI.==13));
   x14=selif(D,(DI.==14));
   x15=selif(D,(DI.==15));
   x16=selif(D,(DI.==16));
   x21=selif(D,(DI.==21));
   x22=selif(D,(DI.==22));
   x23=selif(D,(DI.==23));
   x24=selif(D,(DI.==24));
   x25=selif(D,(DI.==25));
   x26=selif(D,(DI.==26));   
   x31=selif(D,(DI.==31));
   x32=selif(D,(DI.==32));
   x33=selif(D,(DI.==33));
   x34=selif(D,(DI.==34));
   x35=selif(D,(DI.==35));
   x36=selif(D,(DI.==36));
   x41=selif(D,(DI.==41));
   x42=selif(D,(DI.==42));
   x43=selif(D,(DI.==43));
   x44=selif(D,(DI.==44));
   x45=selif(D,(DI.==45));
   x46=selif(D,(DI.==46));
   x51=selif(D,(DI.==51));
   x52=selif(D,(DI.==52));
   x53=selif(D,(DI.==53));
   x54=selif(D,(DI.==54));
   x55=selif(D,(DI.==55));
   x56=selif(D,(DI.==56));
   x61=selif(D,(DI.==61));
   x62=selif(D,(DI.==62));
   x63=selif(D,(DI.==63));
   x64=selif(D,(DI.==64));
   x65=selif(D,(DI.==65));
   x66=selif(D,(DI.==66));
/* joint distribution */
   m=rows(x11)~rows(x12)~rows(x13)~rows(x14)~rows(x15)~rows(x16)|
     rows(x21)~rows(x22)~rows(x23)~rows(x24)~rows(x25)~rows(x26)|
     rows(x31)~rows(x32)~rows(x33)~rows(x34)~rows(x35)~rows(x36)|
     rows(x41)~rows(x42)~rows(x43)~rows(x44)~rows(x45)~rows(x46)|
     rows(x51)~rows(x52)~rows(x53)~rows(x54)~rows(x55)~rows(x56)|
     rows(x61)~rows(x62)~rows(x63)~rows(x64)~rows(x65)~rows(x66);
   joint=m/n;
/* Gibbs sampling */
      a=1; b=1;           /* beginning from (1,1) you could change depending on data. */
      temp=zeros(n,1);
      i=1;
      do while i<=n;
         a=usample(seqa(1,1,rows(joint)),joint[.,b]/sumc(joint[.,b]));
         b=usample(seqa(1,1,cols(joint)),joint[a,.]'/sumc(joint[a,.]'));
         if a==1 and b==1;
            index=ceil(rows(x11)*rndu(1,1));
            temp[i]=x11[index];
         elseif a==1 and b==2;
            index=ceil(rows(x12)*rndu(1,1));
            temp[i]=x12[index];
         elseif a==1 and b==3;
            index=ceil(rows(x13)*rndu(1,1));
            temp[i]=x13[index];         
         elseif a==1 and b==4;
            index=ceil(rows(x14)*rndu(1,1));
            temp[i]=x14[index]; 
         elseif a==1 and b==5;
            index=ceil(rows(x15)*rndu(1,1));
            temp[i]=x15[index];             
         elseif a==1 and b==6;
            index=ceil(rows(x16)*rndu(1,1));
            temp[i]=x16[index];          
         elseif a==2 and b==1;
            index=ceil(rows(x21)*rndu(1,1));
            temp[i]=x21[index];         
         elseif a==2 and b==2;
            index=ceil(rows(x22)*rndu(1,1));
            temp[i]=x22[index];         
         elseif a==2 and b==3;
            index=ceil(rows(x23)*rndu(1,1));
            temp[i]=x23[index];         
         elseif a==2 and b==4;
            index=ceil(rows(x24)*rndu(1,1));
            temp[i]=x24[index];
         elseif a==2 and b==5;
            index=ceil(rows(x25)*rndu(1,1));
            temp[i]=x25[index];         
         elseif a==2 and b==6;
            index=ceil(rows(x26)*rndu(1,1));
            temp[i]=x26[index];     
         elseif a==3 and b==1;
            index=ceil(rows(x31)*rndu(1,1));
            temp[i]=x31[index];      
         elseif a==3 and b==2;
            index=ceil(rows(x32)*rndu(1,1));
            temp[i]=x32[index];        
         elseif a==3 and b==3;
            index=ceil(rows(x33)*rndu(1,1));
            temp[i]=x33[index];    
         elseif a==3 and b==4;
            index=ceil(rows(x34)*rndu(1,1));
            temp[i]=x34[index];
         elseif a==3 and b==5;
            index=ceil(rows(x35)*rndu(1,1));
            temp[i]=x35[index];  
         elseif a==3 and b==6;
            index=ceil(rows(x36)*rndu(1,1));
            temp[i]=x36[index];   
         elseif a==4 and b==1;
            index=ceil(rows(x41)*rndu(1,1));
            temp[i]=x41[index]; 
         elseif a==4 and b==2;
            index=ceil(rows(x42)*rndu(1,1));
            temp[i]=x42[index];          
         elseif a==4 and b==3;
            index=ceil(rows(x43)*rndu(1,1));
            temp[i]=x43[index];          
         elseif a==4 and b==4;
            index=ceil(rows(x44)*rndu(1,1));
            temp[i]=x44[index];
         elseif a==4 and b==5;
            index=ceil(rows(x45)*rndu(1,1));
            temp[i]=x45[index];       
         elseif a==4 and b==6;
            index=ceil(rows(x46)*rndu(1,1));
            temp[i]=x46[index];
         elseif a==5 and b==1;
            index=ceil(rows(x51)*rndu(1,1));
            temp[i]=x51[index]; 
         elseif a==5 and b==2;
            index=ceil(rows(x52)*rndu(1,1));
            temp[i]=x52[index];         
         elseif a==5 and b==3;
            index=ceil(rows(x53)*rndu(1,1));
            temp[i]=x53[index];        
         elseif a==5 and b==4;
            index=ceil(rows(x54)*rndu(1,1));
            temp[i]=x54[index];
         elseif a==5 and b==5;
            index=ceil(rows(x55)*rndu(1,1));
            temp[i]=x55[index];          
         elseif a==5 and b==6;
            index=ceil(rows(x56)*rndu(1,1));
            temp[i]=x56[index];
         elseif a==6 and b==1;
            index=ceil(rows(x61)*rndu(1,1));
            temp[i]=x61[index]; 
         elseif a==6 and b==2;
            index=ceil(rows(x62)*rndu(1,1));
            temp[i]=x62[index];         
         elseif a==6 and b==3;
            index=ceil(rows(x63)*rndu(1,1));
            temp[i]=x63[index];       
         elseif a==6 and b==4;
            index=ceil(rows(x64)*rndu(1,1));
            temp[i]=x64[index];
         elseif a==6 and b==5;
            index=ceil(rows(x65)*rndu(1,1));
            temp[i]=x65[index];        
         else;
            index=ceil(rows(x66)*rndu(1,1));
            temp[i]=x66[index];  
         endif;           
         i=i+1;
      endo;
   r=ones(n,1)+temp;
   S=ST*cumprodc(r);
   retp(ST|S);
endp;

proc usample(x,p);
   local n,u,y;
   n=rows(x);
   if n>10;
      errorlog "ERROR: This function can sample only n<=10.";
      retp(".");
   endif;
   if rows(x)/=rows(p);
      errorlog "ERROR: Rows of x and its p must be the same.";
      retp(".");
   endif;
   _fcmptol=1e-2;
   if fne(sumc(p),1);
      errorlog "ERROR: Sum of the probabilities must be 1.";
      retp(".");
   endif;
   u=rndu(1,1);
   if u>=0 and u<p[1];
      y=x[1];
   elseif u>=p[1] and u<(p[1]+p[2]);
      y=x[2];
   elseif u>=p[2] and u<(p[1]+p[2]+p[3]);
      if rows(x)<3;
         y=x[2];
      else;
         y=x[3];
      endif;
   elseif u>=p[3] and u<(p[1]+p[2]+p[3]+p[4]);
      if rows(x)<4;
         y=x[3];
      else;
         y=x[4];
      endif;
   elseif u>=p[4] and u<(p[1]+p[2]+p[3]+p[4]+p[5]);
      if rows(x)<5;
         y=x[4];
      else;
         y=x[5];
      endif;
   elseif u>=p[2] and u<(p[1]+p[2]+p[3]+p[4]+p[5]+p[6]);
      if rows(x)<6;
         y=x[5];
      else;
         y=x[6];
      endif;
   elseif u>=p[2] and u<(p[1]+p[2]+p[3]+p[4]+p[5]+p[6]+p[7]);
      if rows(x)<7;
         y=x[6];
      else;
         y=x[7];
      endif;
   elseif u>=p[2] and u<(p[1]+p[2]+p[3]+p[4]+p[5]+p[6]+p[7]+p[8]);
      if rows(x)<8;
         y=x[7];
      else;
         y=x[8];
      endif;
   elseif u>=p[2] and u<(p[1]+p[2]+p[3]+p[4]+p[5]+p[6]+p[7]+p[8]+p[9]);
      if rows(x)<9;
         y=x[8];
      else;
         y=x[9];
      endif;
   else;
      y=x[10];
   endif;
   retp(y);
endp;