
/****Copyright (c) 2002  Alexandros V. Gerbessiotis. 
 **** All rights reserved.                               
 ****Permission to use, copy, modify, and distribute this software,
 ****and to incorporate it, in whole or in part, into 
 ****other software,is hereby granted without fee, provided that
 ****  (1) the above copyright notice and this permission 
 ****      notice appear in all copies of the source code, and 
 ****      the above copyright notice appear in clearly visible 
 ****      form on all supporting documentation and distribution 
 ****      media, and
 ****     (2) any redistribution of the software, in original
 ****         or modified form, be without fee and subject to
 ****         these same conditions.
 ****   No guarantee, whatsoever, is offered that the code works as
 ****   advertised or is absent of any, even damaging,
 ****   side-effects. Use this code at your own personal risk.
 ****/
#include "option.h"
#define OPTION_DEBUG
#undef OPTION_DEBUG

double option(double K, double T, double S, double r, int N,  
             double sigma, double dx, double dividend, 
             double *rstock, double *roption, int opt)
{

 double dt,n,p,discount,invexpdx;
 double dxu,dxd,pu,pm,pd,v;
 double *St,*C;
 int i,j,k;
 double delta,gamma;
 double a1,a2,b2,c1,c2;
 double dpd,dpu,dpm,oneoveru;

 St = (double *) malloc((2*N+1)*sizeof(double));
#ifdef SANITY_CHECK
 check_if_null((void *)St,"St","option");
#endif
 C  = (double *) malloc((2*N+1)*sizeof(double));
#ifdef SANITY_CHECK
 check_if_null((void *)St,"C","option");
#endif


 dt = (double) T / (double) N;
 n = (double) N;
 v = r - dividend - 0.5*sigma*sigma;
 invexpdx=1.0/exp(dx);
 pu  = 0.5 * ((sigma*sigma*dt+v*v*dt*dt)/(dx*dx) +v*dt/dx);
 pd  = 0.5 * ((sigma*sigma*dt+v*v*dt*dt)/(dx*dx) -v*dt/dx);
 pm  = 1-(sigma*sigma*dt+v*v*dt*dt)/(dx*dx); 
 discount = exp(-r*dt);

#ifdef OPTION_DEBUG
 fprintf(stdout, "dt=%f n=%f v=%f expdx=%f iexpdx=%f (p)=(%f,%f,%f) discount=%f\n",dt,n,v,1.0/invexpdx,invexpdx,pu,pm,pd,discount);
 fflush(stdout);
#endif


 St[0] = S * exp(N*dx);
 for(j=1;j<=2*N;j++)
   St[j] = St[j-1]*invexpdx;
#ifdef OPTION_DEBUG
 PRINTMATRIX("St","Init",St,(2*N+1));
#endif

 if ((opt == AMERICANCALL) || (opt == EUROPEANCALL)) {
    for(j=0;j<=2*N;j++)
      C[j] = MAX((0.0), (St[j]-K));
 }    
 else {  /* It is an AMERICANPUT or a EUROPEANPUT */
    for(j=0;j<=2*N;j++)
      C[j] = MAX((0.0), (K-St[j]));
 }
#ifdef OPTION_DEBUG
 PRINTMATRIX("C ","Init",C,(2*N+1)); NEWLINE;
#endif

 oneoveru= invexpdx;
 dpd=discount*pd;
 dpu=discount*pu;
 dpm=discount*pm;
if ((opt == EUROPEANCALL) ) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<=2*i;j++) {
      C[j] = dpu*C[j] + dpm*C[j+1]+dpd*C[j+2];
      St[j]=St[j]*oneoveru;
   }
#ifdef OPTION_DEBUG
/*
  PRINTMAT("St",i,St,2*i+1); 
*/
  PRINTMAT("C ",i,C,2*i+1); 
#endif

 }
}
if ((opt == AMERICANPUT ) ) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<=2*i;j++) {
      C[j] = dpu*C[j] + dpm*C[j+1]+dpd*C[j+2];
      St[j]=St[j]*oneoveru;
      C[j] = MAX( C[j], (K-St[j]));
   }
#ifdef OPTION_DEBUG
  PRINTMAT("St",i,St,2*i+1); 
  PRINTMAT("C ",i,C,2*i+1); 
#endif
 }
}
if ((opt == AMERICANCALL) ) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<=2*i;j++) {
      C[j] = dpu*C[j] + dpm*C[j+1]+dpd*C[j+2];
      St[j]=St[j]*oneoveru;
      C[j] = MAX( C[j], (St[j]-K));
   }
 }
}
if ((opt == EUROPEANPUT ) ) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<=2*i;j++) {
      C[j] = dpu*C[j] + dpm*C[j+1]+dpd*C[j+2];
      St[j]=St[j]*oneoveru;
   }
 }
}

*rstock= St[0]; *roption= C[0];
return(*roption);
}



