new; cls;
/* Block Bootstrap Asian Call & Put Options            */
/*     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
;
data=P[2:752];
n=rows(data);
blength=10;
ST=data[n];

K=20;
r=0.03;
t=2.9;
times=5000;


print/lz "      n=" n;
print/lz "blength=" blength;
print/lz "  times=" times;
print "C=" BBasianC(ST,data,blength,K,r,t,times);
print "P=" BBasianP(ST,data,blength,K,r,t,times);
print;
print "Another version:";
print "C=" BBasianC2(ST,data,blength,r,t,times);
print "P=" BBasianP2(ST,data,blength,r,t,times);
print;
print "Geometric Mean Asian:";
print "C=" BBgasianC(ST,data,blength,K,r,t,times);
print "P=" BBgasianP(ST,data,blength,K,r,t,times);
print;
print "Another version of Geometric Mean:";
print "C=" BBgasianC2(ST,data,blength,r,t,times);
print "P=" BBgasianP2(ST,data,blength,r,t,times);



/*
**  bbasian.txt - Block Bootstrap Asian Call & Put Options.
** (C) Copyright 2004 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Asian option prices by block-bootstrapped simulation
**              from real past data.
**
**  Format:     C=BBasianC(ST,data,blength,K,r,t,times);
**              P=BBasianP(ST,data,blength,K,r,t,times);
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**             
**              blength  scalar,  an integer of block length (a divisor of n-1)
**
**              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
**
*/

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

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

/*
**  bbasian.txt - Another version of Block Bootstrap Asian Call & Put Options.
** (C) Copyright 2004 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Asian option prices by block-bootstrapped simulation
**              from real past data.
**
**  Format:     C=BBasianC2(ST,data,blength,r,t,times);
**              P=BBasianP2(ST,data,blength,r,t,times);
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**             
**              blength  scalar,  an integer of block length (a divisor of n-1)
**
**              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
**
*/

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

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

/*
**  bbasian.txt - Block Bootstrap Geometric Asian Call & Put Options.
** (C) Copyright 2004 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Geometric Asian option prices by block-bootstrapped simulation
**              from real past data.
**
**  Format:     C=BBgasianC(ST,data,blength,K,r,t,times);
**              P=BBgasianP(ST,data,blength,K,r,t,times);
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**             
**              blength  scalar,  an integer of block length (a divisor of n-1)
**
**              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
**
*/

proc BBgasianC(ST,data,blength,K,r,t,times);
   local n,AV,i,S,C;
   n=rows(data);
   AV=zeros(times,1);
   i=1;
   do while i<=times;
      S=BBsampler2(ST,data,blength);
      AV[i]=geomean(S[2:n]);
      i=i+1;
   endo;
   C=exp(-r*t)*meanc(maxc(((AV-K)~zeros(times,1))'));   /* We use arithmetic one here. */
   retp(C);
endp;

proc BBgasianP(ST,data,blength,K,r,t,times);
   local n,AV,i,S,P;
   n=rows(data);
   AV=zeros(times,1);
   i=1;
   do while i<=times;
      S=BBsampler2(ST,data,blength);
      AV[i]=geomean(S[2:n]);
      i=i+1;
   endo;
   P=exp(-r*t)*meanc(maxc(((K-AV)~zeros(times,1))'));   /* we use arithmetic one here. */
   retp(P);
endp;

/*
**  bbasian.txt - Another version of Block Bootstrap Geometric Asian Call & Put Options.
** (C) Copyright 2004 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Geometric Asian option prices by block-bootstrapped simulation
**              from real past data.
**
**  Format:     C=BBgasianC2(ST,data,blength,r,t,times);
**              P=BBgasianP2(ST,data,blength,r,t,times);
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**             
**              blength  scalar,  an integer of block length (a divisor of n-1)
**
**              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
**
*/

proc BBgasianC2(ST,data,blength,r,t,times);
   local n,V,i,S,C;
   n=rows(data);
   V=zeros(times,1);
   i=1;
   do while i<=times;
      S=BBsampler2(ST,data,blength);
      V[i]=S[n]-geomean(S[2:n]);
      i=i+1;
   endo;
   C=exp(-r*t)*meanc(maxc((V~zeros(times,1))'));
   retp(C);
endp;

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

/*
**  bbsampler.txt - Time Series Change Amount Block Sampler from real data.
** (C) Copyright 2004 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Samples just one block-bootstrap representative data from real data.
**
**  Format:     S=BBsampler(ST,data,blength);
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**             
**              blength  scalar,  an integer of block length (a divisor of n-1)
**
**
**  Output:     S        vector,  n x 1 vector of block-bootstrapped representative data
**
**
**  Notice:     Some of S could go below zero.
**
*/

proc BBsampler(ST,data,blength);
   local n,delS,col,m,index,del0t,S;
   n=rows(data);
   delS=data[2:n]-data[1:n-1];
   col=(n-1)/blength;
   if round(col)-col/=0;
      errorlog "ERROR: Parameter blength must be a divisor of n-1.";
      retp(".");
   endif;  
   m=reshape(delS,col,blength)';
   index=ceil(col*rndu(col,1));
   del0t=m[.,index];
   S=ST+cumsumc(reshape(del0t',n-1,1));
   retp(ST|S);
endp;


/*
**  bbsampler.txt - Time Series Change Rate Block Sampler from real data.
** (C) Copyright 2004 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Samples just one block-bootstrap representative data from real data.
**
**  Format:     S=BBsampler2(ST,data,blength);
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**             
**              blength  scalar,  an integer of block length (a divisor of n-1)
**
**
**  Output:     S        vector,  n x 1 vector of block-bootstrapped representative data
**
*/

proc BBsampler2(ST,data,blength);
   local n,d,col,m,index,d0t,S;
   n=rows(data);
   d=ones(n-1,1)+(ln(data[2:n])-ln(data[1:n-1]));
   col=(n-1)/blength;
   if round(col)-col/=0;
      errorlog "ERROR: Parameter blength must be a divisor of n-1.";
      retp(".");
   endif;  
   m=reshape(d,col,blength)';
   index=ceil(col*rndu(col,1));
   d0t=m[.,index];
   S=ST*cumprodc(reshape(d0t',n-1,1));
   retp(ST|S);
endp;

proc geomean(x);
   local n;
   n=rows(x);
   retp(exp(1/n*sumc(ln(x))));
endp;