/*
** Version 01
** Date : 20/07/2001                                        
** --------------------------------------------------------------------------
** Copyright K.Cuthbertson and D. Nitzsche
** "Financial Engineering:Derivatives and Risk Manangement" - J. Wiley 2001
** 
**   To calculate the Black-Scholes 'Greeks' for European option which pays 
**   'dividends' and to graph the 'Greeks' against S.  
**   The formulae used for the Greeks can be found in Appendix Cht 9 and 
**   Fig 9.2 and 9.4 plus some other graphs are produced below. 
**
*/

new ; cls ; format /m1/rd 12,6;  output on ; screen on ;

/*
** ofile = "c:\\docs\\arun\\out.out" ;
** output file = ^ofile reset ; 
*/

/* --------------------------------------------------------------------------------
** USER INPUT: 
** USER MAY CHANGE FIGURES IN THE NEXT SECTION 
** -------------------------------------------------------------------------------- */
k     = 50 ;        @ strike price @
s0    = 50 ;        @ price underlying asset @   
r     = 10 ;        @ interest rate, contin. comp, percent p.a. @
d     = 0.0 ;       @ 'dividend yield', contin. comp, percent p.a. @
sigma = 20  ;       @ standard deviation, annual, percent p.a. @
tau   = 0.5 ;       @ time to maturity, years @

_type_option = "call";        @ ---- User chooses "call" or "put" in 'lower case'- see procedures below ---- @
 
/* ------------------------------------------  END OF USER INPUT ----------------------------------------- */

sigma  = sigma./100  ;
r      = r./100      ;		
d      = d./100     ;

" ---------------------------------------- Data and Outputs   ------------------------------------------  " ;
?;
"   Value of underlying asset (s)                     "  s0    ;  
"   Strike price (k)                                  "  k     ; 
"   Time to maturity (tau, years or fraction of year) "  tau   ; 
"   Volatility of underlying (sigma, propn. )         "  sigma ; 
"   Risk free rate (r, contin. comp, propn)           "  r     ;
"   Dividend yield (d, contin. comp, propn)           "  d     ; 
?;
"   Initial value of the call and put premia          "  bscall(k,s0,sigma,tau,r,d)~bsput(k,s0,sigma,tau,r,d); 
"   Initial call delta ($ per $)                      "  bs_delta(k,s0,sigma,tau,r,d);
"   Initial call theta ($ per day )                   "  bs_theta(k,s0,sigma,tau,r,d)./365 ;
"   Initial call vega ($ per 1% change)               "  bs_vega (k,s0,sigma,tau,r,d)/100 ;
"   Initial call gamma ($, [per $]^2 )                "  bs_gamma(k,s0,sigma,tau,r,d);
"   Initial call rho ($ per 1% change)                "  bs_rho  (k,s0,sigma,tau,r,d)./100 ;
" -------------------------------------------------------------------------------------------------------  " ;

_type_option = "put";        @ User chooses "call" or "put" in 'lower case'- see procedures below @

?;
"   Initial value of the call and put premia          "  bscall(k,s0,sigma,tau,r,d)~bsput(k,s0,sigma,tau,r,d); 
"   Initial put delta ($ per $)                       "  bs_delta(k,s0,sigma,tau,r,d);
"   Initial put theta ($ per day )                    "  bs_theta(k,s0,sigma,tau,r,d)./365 ;
"   Initial put vega ($ per 1% change)                "  bs_vega (k,s0,sigma,tau,r,d)./100 ;
"   Initial put gamma ($, [per $]^2 )                 "  bs_gamma(k,s0,sigma,tau,r,d);
"   Initial put rho ($ per 1% change)                 "  bs_rho  (k,s0,sigma,tau,r,d)./100;
" -------------------------------------------------------------------------------------------------------  " ;


@ ----------------------  Setting up data for graphs of Greeks for a CALL ------------------------------ @ 

_type_option = "call";        @ User chooses "call" or "put" in 'lower case'- see procedures below @


s      = seqa(s0-20,1,40);    @ Sequence 40 values for s, starting at s0-20 and ending at s0+20 @ 
xx     = zeros(40,5);         @ The 5 'Greeks' are to be stored to in vector xx @ 

@ ------------------------------ Activate the procedures for the Greeks ----------------------- @
xx[.,1] = bs_delta(k,s,sigma,tau,r,d);
xx[.,2] = bs_theta(k,s,sigma,tau,r,d);
xx[.,3] = bs_vega (k,s,sigma,tau,r,d);
xx[.,4] = bs_gamma(k,s,sigma,tau,r,d);
xx[.,5] = bs_rho  (k,s,sigma,tau,r,d);


@ --------------------------------  Set up the multiple graphics output ------------------------ @
library pgraph;
graphset;

begwind;
window(2,2,0);
  axmargin(1,1,1,1);
  
  xlabel("Price of underlying");
  ylabel("Delta of Call");
  xy(s,xx[.,1]);
  
nextwind;
  ylabel("Vega of Call");
  xy(s,xx[.,3]);
  
nextwind;
  ylabel("Gamma of Call");
  xy(s,xx[.,4]);
  
nextwind;
  ylabel("Rho of Call");
  xy(s,xx[.,5]);

endwind ;

@ ---------------------------------- Separate graph for theta ---------------------------------  @

 ylabel("Theta of Call");
  xy(s,xx[.,2]);

@ ----------------------------- Graph delta and Gamma together --------------------------------- @
z = xx[.,1]~10.*xx[.,4] ;
 ylabel("Delta and Gamma(x10) of Call");
  xy(s,z);

@ ---------------  Setting up data for 3-D graph of gamma against s and tau for a put ----------- @ 

s1   = seqa(40,.1,200);
tau1 = seqa(tau,-1/365,130);   @ tau,initially= 0.5. New tau1 is from 0.5 to '0.5-130/365' >0  @ 

_type_option = "put" ;         @ Switch to gamma for a put option @

xxx = zeros(200,130);

i=1;
do until i>130;
  xxx[.,i]=bs_gamma(k,s1,sigma,tau1[i],r,d);
  i=i+1;
endo;

@ ----------------------  3-D Graph for the Gamma v/s Underlying and Maturity --------------------- @

 graphset; 
  xlabel("Maturity");
  ylabel("Price of underlying");
  zlabel("Gamma of put option");
     _pdate="";
     _pframe=0;
surface(10.*tau1',s1,xxx);

@ -----------------------------------------  End of graphics output --------------------------------- @


@ ----------------------------------------------   Procedures --------------------------------------- @

@ ------------------------ Black Scholes (with cont. comp.dividends) see appendix 9.1 --------------- @
  
proc bscall(k,s,sigma,tau,r,d);
  local d1,d2,c;
  
  d1 = ( ln( s./k ) + ( r - d + sigma^2./2).*tau )./ ( sigma.*sqrt(tau) );
  d2 = d1 - sigma.*sqrt(tau);
  
  c=s.*exp(-d.*tau).*cdfn(d1) -k.*exp(-r.*tau).*cdfn(d2);
  
  retp(c);
endp;


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

@ ------------------------------------ The Greeks, see appendix 9.1 --------------------------- @

@ -------------------------------------------------- Delta  ----------------------------------- @

proc bs_delta(k,s,sigma,tau,r,d);
  local d1,delta;
  
   d1 = ( ln( s./k ) + ( r - d + sigma^2./2).*tau )./ ( sigma.*sqrt(tau) );
   
        if (lower(_type_option) $== "call");      delta = cdfn(d1);
    elseif (lower(_type_option) $== "put");       delta = cdfn(d1)-1;
  
  else;
    errorlog "Error : _type_option = call or _type_option = put";
    retp(error(0));
  endif;
  
  retp(delta);
endp;


@ ------------------------------------------ Gamma ----------------------------------------------- @
@ -------------- Note that 'pdfn(d1)' is the Gauss command for N'(d1) in the text ---------------- @

proc bs_gamma(k,s,sigma,tau,r,d);

  local d1,gam;                     @ Cannot use 'gamma', as this is Gauss internal command @
  
  d1 = ( ln( s./k ) + ( r - d + sigma^2./2).*tau )./ ( sigma.*sqrt(tau) ); 
  
 
        if (lower(_type_option) $== "call");      gam = pdfn(d1)./(s.*sigma.*sqrt(tau));
    elseif (lower(_type_option) $== "put" );      gam = pdfn(d1)./(s.*sigma.*sqrt(tau));
          
  else;
    errorlog " Error : _type_option = call or _type_option = put " ;
    retp(error(0));
  endif;
  
  retp(gam);
endp;

@ ------------------------------------------- Theta -------------------------------------------- @

proc bs_theta(k,s,sigma,tau,r,d);
  local d1,d2,theta;
  
  d1 = ( ln( s./k ) + ( r - d + sigma^2./2).*tau )./ ( sigma.*sqrt(tau) );
  d2 = d1 - sigma.*sqrt(tau);
   
  if (lower(_type_option) $== "call");
  
  theta = -( s.*sigma.*pdfn(d1).*exp(-d.*tau) )./(2*sqrt(tau)) - k.*r.*exp(-r.*tau).*cdfn(d2)
             + d.*s.*cdfn(d1).*exp(-d.*tau);
            
  elseif (lower(_type_option) $== "put");
   
  theta = -( s.*sigma.*pdfn(d1).*exp(-d.*tau) )./(2*sqrt(tau)) + k.*r.*exp(-r.*tau).*cdfn(-d2)
             + d.*s.*cdfn(-d1).*exp(-d.*tau);
      
  else;
    errorlog "Error : _type_option = call or _type_option = put";
    retp(error(0));
  endif;
  
  retp(theta);
  
endp;

@ ------------------------------------------- Vega ------------------------------------------------ @

/* Vega */

proc bs_vega(k,s,sigma,tau,r,d);
  local d1,vega ;
  
  d1 = ( ln( s./k ) + ( r - d + sigma^2./2).*tau )./ ( sigma.*sqrt(tau) );
  
 /* d1=(ln(s./(k.*exp(-r.*tau))))./(sigma.*sqrt(tau))+0.5*sigma.*sqrt(tau); */
  
      if (lower(_type_option) $== "call");      vega = s.*sqrt(tau).*pdfn(d1);
  elseif (lower(_type_option) $== "put" );      vega = s.*sqrt(tau).*pdfn(d1);
  
  else;
    errorlog "Error : _type_option = call or _type_option = put";
    retp(error(0));
  endif;
  
  retp(vega);
  
endp;

@ ----------------------------------------------- Rho --------------------------------------------- @

/* Rho */

proc bs_rho(k,s,sigma,tau,r,d);
  local d1,d2,rho;
  
  d1 = ( ln( s./k ) + ( r - d + sigma^2./2).*tau )./ ( sigma.*sqrt(tau) );
  d2 = d1 - sigma.*sqrt(tau);
  
  /* d2=(ln(s./(k.*exp(-r.*tau))))./(sigma.*sqrt(tau))-0.5*sigma.*sqrt(tau); */
  
      if (lower(_type_option) $== "call");          rho=tau.*k.*exp(-r.*tau).*cdfn(d2);
  elseif (lower(_type_option) $== "put" );          rho=tau.*k.*exp(-r.*tau).*(cdfn(d2)-1);
 
  else;
    errorlog "Error : _type_option = call or _type_option = put";
    retp(error(0));
  endif;
  
  retp(rho);
endp;


@ --------------------------------- End of procedures -------------------------------------- @

end ;