new; cls;
/* Empirical Monte Carlo alpha-percentile Call Options */
/* restricted by Scenarios (using first 3 scenarios)   */
/*     Data from Sheldon M. Ross(1999) Table 10.5 in   */
/*     "An Introduction to Mathematical Finance"       */
/*                             (Wait for a 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];
nn=252;
K=20;
r=0.03;
t=1;
alpha=0.75;
times=200;                                 /* Try 5000 */



print/lz "  times=" times;
print/lz "  alpha=" alpha;
print "floating strike:";
print "C=" SEMCalphapC(ST,data,nn,r,t,alpha,times);
print;
print "fixed strike:";
print "C=" SEMCalphapCf(ST,data,nn,K,r,t,alpha,times);




/*
**  semcalphap.txt - Restriction Scenarios Emprical Monte Carlo alpha-percentile(Floating Strike) Call Options.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates alpha-percentile option prices by Emprical Monte Carlo simulation
**              restricted by Scenarios from real past data.
**
**  Format:     C=SEMCalphapC(ST,data,nn,r,t,alpha,times);
**              
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**             
**              nn       scalar,  number of business times(days,weeks,etc..) until maturity
**
**              r        scalar,  risk-free interest rate
**
**              t        scalar,  maturity (the same as data period)
**
**              alpha    scalar,  percentile ( 0 <= alpha <= 1 )
**
**              times    scalar,  times of simulation
**
**
**  Output:     C        scalar,  call option price
**
**
**  Notice:     It takes very long time to run. Try some small 'times' first.
**
*/

proc SEMCalphapC(ST,data,nn,r,t,alpha,times);
   local OV,i,S,STT,C;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=scenarioemcsampler(ST,data,nn);
      STT=S[nn];
      S=S[2:nn];
      S=sortc(S,1);
      if round((nn-1)*alpha)==0;
          OV[i]=STT-S[1];
      else;
          OV[i]=STT-S[round((nn-1)*alpha)];
      endif;
      i=i+1;
   endo;
   C=exp(-r*t)*meanc(maxc((OV~zeros(times,1))'));
   retp(C);
endp;

/*
**  semcalphap.txt - Restriction Scenarios Emprical Monte Carlo alpha-percentile(Fixed Strike) Call Options.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates alpha-percentile option prices by Emprical Monte Carlo simulation
**              restricted by Scenarios from real past data.
**
**  Format:     C=SEMCalphapCf(ST,data,nn,r,K,t,alpha,times);
**              
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**             
**              nn       scalar,  number of business times(days,weeks,etc..) until maturity
**
**              K        scalar,  strike price
**
**              r        scalar,  risk-free interest rate
**
**              t        scalar,  maturity (the same as data period)
**
**              alpha    scalar,  percentile ( 0 <= alpha <= 1 )
**
**              times    scalar,  times of simulation
**
**
**  Output:     C        scalar,  call option price
**
**
**  Notice:     It takes very long time to run. Try some small 'times' first.
**
*/

proc SEMCalphapCf(ST,data,nn,K,r,t,alpha,times);
   local OV,i,S,C;
   OV=zeros(times,1);
   i=1;
   do while i<=times;
      S=scenarioemcsampler(ST,data,nn);
      S=S[2:nn];
      S=sortc(S,1);
      if round((nn-1)*alpha)==0;
          OV[i]=S[1]-K;
      else;
          OV[i]=S[round((nn-1)*alpha)]-K;
      endif;
      i=i+1;
   endo;
   C=exp(-r*t)*meanc(maxc((OV~zeros(times,1))'));
   retp(C);
endp;


/* This is an example sampler. */

proc scenarioemcsampler(ST,data,nn); 
   local S;
   S=emcsampler2(ST,data,nn);
   do while (scenario01(S,10,10) or scenario02(S,2,0.5) or scenario03(S,0.10,20));
      S=emcsampler2(ST,data,nn);
   endo;
   retp(S);
endp;


/*
**  scenario01.txt - restriction scenario(set max consecutive up and down days).
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Get the status whether or not path S is resticted.
**
**  Format:     status=scenario01(S,maxupdays,maxdowndays);
**
**  Input:      S              vector,  a path to be determined whether or not it is used
**
**              maxupdays      scalar,  maximum consecutive upward(rising) days
**
**              maxdowndays    scalar,  maximum consecutive downward(falling) days
**             
**
**  Output:     status         scalar,  1 if restricted, 0 otherwise
**
*/
proc scenario01(S,maxupdays,maxdowndays);
   local n,delta,signs,i,count,result;
   n=rows(S);
   delta=S[2:n]-S[1:n-1];
/* looking upward */
   signs=(delta.>0);
   result=zeros(n-1,1);
   count=0;
   i=1;
   do while i<=n-1;
      if signs[i]==0;
         count=0;
         result[i]=count;
      else;
         count=count+1;
         result[i]=count;
      endif;
      i=i+1;
   endo;
   if maxc(result)>maxupdays;
      retp(1);
   endif;
/* looking downward */
   signs=(delta.<0);
   result=zeros(n-1,1);
   count=0;
   i=1;
   do while i<=n-1;
      if signs[i]==0;
         count=0;
         result[i]=count;
      else;
         count=count+1;
         result[i]=count;
      endif;
      i=i+1;
   endo;
   if maxc(result)>maxdowndays;
      retp(1);
   endif;
   retp(0);
endp;

/*
**  scenario02.txt - restriction scenario(set maximum amplitude ratios from the initial price).
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Get the status whether or not path S is resticted.
**
**  Format:     status=scenario02(S,maxupratio,maxdownratio);
**
**  Input:      S              vector,  a path to be determined whether or not it is used
**
**              maxupratio     scalar,  maximum upswing amplitude ratio from S0
**
**              maxdownratio   scalar,  maximum downswing amplitude rate from S0
**             
**
**  Output:     status         scalar,  1 if restricted, 0 otherwise
**
**  Example     if maxupratio=2 and maxdownratio=0.5, it means 0.5*S0 <= S <= 2*S0 leading to 0.
*/
proc scenario02(S,maxupratio,maxdownratio);
/* looking upswing */
   if maxc(S)>maxupratio*S[1];
      retp(1);
   endif;
/* looking downswing */
   if minc(S)<maxdownratio*S[1];
      retp(1);
   endif;
   retp(0);
endp;

/*
**  scenario03.txt - restriction scenario(set maximum drop rate within the first days)).
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Get the status whether or not path S is restricted.
**
**  Format:     status=scenario03(S,maxfirstdroprate,firstdays);
**
**  Input:      S                 vector,  a path to be determined whether or not it is used
**
**              maxfirstdroprate  scalar,  maximum drop rate from S0
**
**              firstdays         scalar,  looking the 'firstdays' days
**             
**
**  Output:     status            scalar,  1 if restricted, 0 otherwise
**
**  Example     if (1-maxfirstdroprate)*S0>S, then 1.
*/
proc scenario03(S,maxfirstdroprate,firstdays);
   if minc(S[1:firstdays+1])<(1-maxfirstdroprate)*S[1];
      retp(1);
   else;
      retp(0);
   endif;
endp;


/*
**  emcsampler.txt - Empirical Monte Carlo Sampler from real data.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Samples just one empirical Monte Carlo representative data from real data.
**
**  Format:     S=emcsampler(ST,data,nn);
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**
**              nn       scalar,  number of resulting values
**             
**
**  Output:     S        vector,  nn x 1 vector of resulting values
**
*/
proc emcsampler(ST,data,nn);
   local x,n,p,y,r,S,j,index;
   x=rndu(nn-1,1);
   data=ln(data[2:rows(data)])-ln(data[1:rows(data)-1]);
   data=sortc(data,1);
   n=rows(data);
   p=seqa(0,1,n)/(n-1);
   y=zeros(nn-1,1);
   j=1;
   do while j<=nn-1;
      if x[j]==1;
         y[j]=maxc(data);
      else;
         index=sumc(x[j].>=p);
         y[j]=data[index]+(data[index+1]-data[index])*rndu(1,1);
      endif;
      j=j+1;
   endo;
   r=ones(nn-1,1)+y;
   S=ST*cumprodc(r);
   retp(ST|S);
endp;

/*
**  emcsampler.txt - Empirical Monte Carlo Sampler from real data.
**                   (Stratified Sampling(LHS) version)
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Samples just one empirical Monte Carlo representative data from real data.
**
**  Format:     S=emcsampler2(ST,data,nn);
**
**  Input:      ST       scalar,  initial value(usually ST from real data)
**
**              data     vector,  n x 1 vector of real past data
**
**              nn       scalar,  number of resulting values
**             
**
**  Output:     S        vector,  nn x 1 vector of resulting values
**
*/
proc emcsampler2(ST,data,nn);
   local d,i,U,x,n,p,y,r,S,j,index;
/* Stratified Sampling(LHS) */
   /* divide between 0 and 1 into n segments and sample from each segment */
   d=1/nn;
   U=zeros(nn,1);
   i=1;
   do while i<=nn;
      U[i]=(d*i-d*(i-1))*rndu(1,1)+d*(i-1);
      i=i+1;
   endo;
   /* check if the first element is zero or not */
   do while U[1]==0;
      U[1]=d*rndu(1,1);
   endo;
   /* randomize the order */
   x=zeros(nn,1);
   i=1;
   do while i<=nn-1;
      index=ceil((nn-(i-1))*rndu(1,1));
      x[i]=U[index];
      if index==1;
         U=U[2:rows(U)];
      elseif index==rows(U);
         U=U[1:rows(U)-1];
      else;
         U=U[1:(index-1) (index+1):rows(U)];
      endif;
      i=i+1;
   endo;
   x[nn]=U;
/* main */
   data=ln(data[2:rows(data)])-ln(data[1:rows(data)-1]);
   data=sortc(data,1);
   n=rows(data);
   p=seqa(0,1,n)/(n-1);
   y=zeros(nn-1,1);
   j=1;
   do while j<=nn-1;
      if x[j]==1;
         y[j]=maxc(data);
      else;
         index=sumc(x[j].>=p);
         y[j]=data[index]+(data[index+1]-data[index])*rndu(1,1);
      endif;
      j=j+1;
   endo;
   r=ones(nn-1,1)+y;
   S=ST*cumprodc(r);
   retp(ST|S);
endp;