new; cls;
/* Geometric Brownian Motion restricted by Scenarios */
S0=100;
nu=0.15;
sig=0.40;
delt=1/252;
nn=252;

S=sGBM(S0,nu,sig,delt,nn);
library pgraph;
graphset;
title("A Path of GBM restricted by Scenarios");
xy(seqa(0,delt,nn+1),S);




proc sGBM(S0,nu,sig,delt,nn); 
   local S;
   S=GBM(S0,nu,sig,delt,nn);
   do while (scenario01(S,10,10) or scenario02(S,2,0.5) or scenario03(S,0.10,20));
      S=GBM(S0,nu,sig,delt,nn);
   endo;
   retp(S);
endp;

proc GBM(S0,nu,sig,delt,nn);
   local e,t,St,i;
   e=rndn(nn,1);            
   t=seqa(0,delt,nn+1);
   St=zeros(nn+1,1);
   St[1]=S0;
   i=1;
   do while i<=nn;
      St[i+1]=exp(nu*delt+sig*e[i]*sqrt(delt))*St[i];
      i=i+1;
   endo;
   retp(St);
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 restricted.
**
**  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 resticted.
**
**  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;