/* 
** Version 01
** Date : 20/07/2001                                        
** --------------------------------------------------------------------------
** Copyright K.Cuthbertson and D. Nitzsche
** "Financial Engineering:Derivatives and Risk Manangement" - J. Wiley 2001
**
**   This prog. calculates the call premium using Black Scholes and the
**   (Cox Ross Rubinstein, CRR ) BOPM model.
**   see Financial Engineering p.229 and Fig 8.9 
** 
**   We calculate a binominal call price for different n's and compare the 
**   outcome with the fixed value given by Black-Scholes       
**                                    
**   In the BOPM   
**        rp   = periodic discrete rate, r=continuously compounded annual rate 
**     Hence R = (1+rp) = exp(r^dt)   
**     and   q = (exp(r^dt)-d)/(u - d)            
**
*/


new ; cls ;
format /m1/rd 12,4;  output on ;

declare external _ucrr ;       @ up value @       @ globals in procs for CRR, to initialize @
declare external _pcrr ;       @ r.n. proby @
declare external _dcrr ;       @ down value @

 " ----------------------------------------------------------------------------------------- " ;
 "                            FILE : Cox Ross Rubinstein BOPM                                " ;
 " ----------------------------------------------------------------------------------------- " ;
 
/*
** " ------------------------------------------------- " ;
** "         Input the variables                       " ;
** "         -------------------                       " ;
** "                                                   " ;
** "   k      = strike price                           " ;   
** "   s      = stock price                            " ;   
** "   sigma  = standard deviation(proportion)         " ;   
** "   tau    = time to maturity(years)                " ;   
** "   r      = contin. compounded rate(p.a.)          " ;  
** "   dt     = time increment                         " ;
** "   nob    = tau/dt = max number of periods in BOPM " ;
** " ------------------------------------------------- " ;
*/
 
@ --------------  USER INPUT: See Financial Engineering p.229 and Fig 8.9  --------------- @ 
s     = 99 ;         
k     = 100 ;                
sigma = 0.20 ;             
r     = 0.10 ;   @ contin. compounded, annual, propn @
tau   = 0.30 ;  

dt    = 0.01 ; 

@ ------------------------------------  END OF USER INPUT ------------------------------------- @

nob   = trunc(tau/dt) ;              @ converts to integer @ 

@ -------------------------- Setting vectors for CALL price using the BS and BOPM ------------- @

call_bs      = zeros(nob,1) ;
call_bopm    = zeros(nob,1) ;
xx           = zeros(nob,1) ;        @ counter for x-axis @

i=1 ;
  do until i > nob ;
     {c_bopm} = coxrosscall( k,s,sigma,tau,r,i ) ;   @ -- Activate the procedures - see below --@
       {c_bs} = bscall( k,s,sigma,tau,r ) ;

     call_bopm[i,1] = c_bopm ; 
     call_bs[i,1]   = c_bs ;  
     xx[i,1] = i ;
  i = i + 1 ;
endo ;


@ ------------------------------------ Print the output --------------------------------------- @  

?; " up=U , risk neutral probability, q and down=D                  "    _ucrr~_pcrr~_dcrr ; 

 " ------------------------------------------------------------------------------------------- " ;
 "   k      = strike price                           "    k     ;       
 "   s      = stock price                            "    s     ;    
 "   sigma  = standard deviation(proportion)         "    sigma ;   
 "   tau    = time to maturity(years)                "    tau   ;  
 "   r      = contin. compounded rate(p.a.)          "    r     ; 
 "   dt     = time increment                         "    dt    ; 
 "   nob    = tau/dt = max number of periods in BOPM "    nob   ;  
 " ------------------------------------------------------------------------------------------ " ;

?; "   Values of BOPM as the number of time steps is increased  " ;
?; "   i /                BS      /         BOPM                " ;
?; " ---------------------------------------------------------  " ; xx~call_bs~call_bopm ;

 

@ ----------------------- Graphical representation BS and BOPM prices -------------------- @

library gauss pgraph ;
graphset ;
screen on ;
  title( "Call Premium: Black Scholes and BOPM Model as n increases " ) ;
     xy( xx , call_bopm~call_bs ) ;


@ ----------------------------- END OF PROGRAM: Procedures follow ------------------------- @         
 
@ ----------------------- Black-Scholes Call and Put ( (without dividends) ---------------- @   

proc bscall(k,s,sigma,tau,r);
     local d1,d2,c ;
        d1=(  ln(s./k) + ( r + sigma^2/2).*tau  )./(  sigma.*sqrt(tau)  ) ;
        d2=d1 - sigma.*sqrt(tau);
     c=s.*cdfn(d1)-k.*exp(-r.*tau).*cdfn(d2);
     retp(c);
endp;

proc bsput(k,s,sigma,tau,r);
     local d1,d2,p ;
        d1=(  ln(s./k) + ( r + sigma^2/2).*tau  )./(  sigma.*sqrt(tau)  ) ;
        d2=d1 - sigma.*sqrt(tau);
     p= k.*exp(-r.*tau).*cdfn(-d2)  -  s.*cdfn(-d1);
     retp(p);
endp;


@ ------------------------   Procs to calculate option premia using CRR ------------------  @ 
@ ----------------------   Call:BOPM for call (r is contin compounded) -------------------  @ 

proc coxrosscall(k,s,sigma,tau,r,n);
  local h,u,d,i,p,q,j,binom,m,aux,sj,c;
  h=maxc(rows(k)|rows(s)|rows(sigma)|rows(tau)|rows(r));
  u=exp(sigma.*sqrt(tau./n));
  d=1./u;
  i=exp(r.*tau./n);
  p=(i-d)./(u-d);
  q=1-p;
  j=seqa(0,1,n+1);
  binom=(n!)./((j!).*((n-j)!));
  sj={};
  j=0;
  do until j>n;
     m=s.*(u^j).*(d^(n-j))-k;
     m=m~zeros(h,1);
     m=maxc(m');
     aux=m.*binom[j+1].*(p^j).*(q^(n-j));
     sj=sj~aux;
     j=j+1;
  endo;
  c=exp(-r.*tau).*sumc(sj');
  _ucrr=u;                           @ --------- globals to initialize IN MAIN PROG ---------- @
  _pcrr=p;
  _dcrr=d;
  retp(c);
endp;

proc coxrossput(k,s,sigma,tau,r,n);
  local h,u,d,i,p,q,j,binom,m,aux,sj,c;
  h=maxc(rows(k)|rows(s)|rows(sigma)|rows(tau)|rows(r));
  u=exp(sigma.*sqrt(tau./n));
  d=1./u;
  i=exp(r.*tau./n);
  p=(i-d)./(u-d);
  q=1-p;
  j=seqa(0,1,n+1);
  binom=(n!)./((j!).*((n-j)!));
  sj={};
  j=0;
  do until j>n;
     m=k-s.*(u^j).*(d^(n-j));
     m=m~zeros(h,1);
     m=maxc(m');
     aux=m.*binom[j+1].*(p^j).*(q^(n-j));
     sj=sj~aux;
     j=j+1;
  endo;
  c=exp(-r.*tau).*sumc(sj');
  _ucrr=u;                             
  _ucrr=u;                          @ --------- globals to initialize IN MAIN PROG ---------- @
  _pcrr=p;
  _dcrr=d;
  retp(c);
endp;

end ;