/****   Copyright (c) 2002  Alexandros V. Gerbessiotis
 ****
 ****   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"


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

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

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


 dt = (double) T / (double) N;
 n = (double) N;
 v = r - 0.5*sigma*sigma;
 dxu = sqrt(sigma*sigma*dt+v*dt*v*dt);
 dxd = -dxu;
 pu  = 0.5 + 0.5 * (v*dt /dxu);
 pd  = 1-pu;
 discount = exp(-r*dt);

#ifdef SANITY_CHECK
 fprintf(stdout, "K=%f T=%d S=%f r=%f N=%d, u=%f d=%f\ndt=%f\
                  p=%f discount=%f\n",K,T,S,r,N,u,d,dt,p,discount);
 fflush(stdout);
#endif


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

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

 oneoveru= 1.0/exp(dxu);
 dpd=discount*pd;
 dpu=discount*pu;
if ((opt == AMERICANCALL) ) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<=i;j++) {
      C[j] = dpu*C[j] + dpd*C[j+1];
      St[j]=St[j]*oneoveru;
      C[j] = MAX( C[j], (St[j]-K));
   }
#ifdef AMERICANOPTION_DEBUG
 PRINTMATRIX("C ","I   ",C,(i+1)); 
#endif
 }
}
if ((opt == AMERICANPUT ) ) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<=i;j++) {
      C[j] = dpu*C[j] + dpd*C[j+1];
      St[j]=St[j]*oneoveru;
      C[j] = MAX( C[j], (K-St[j]));
   }
#ifdef AMERICANOPTION_DEBUG
 PRINTMATRIX("C ","I   ",C,(i+1));
#endif
 }
}
if ((opt == EUROPEANCALL) ) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<=i;j++) {
      C[j] = dpu*C[j] + dpd*C[j+1];
      St[j]=St[j]*oneoveru;
   }
#ifdef EUROPEANOPTION_DEBUG
 PRINTMATRIX("C ","I   ",C,(i+1));
#endif
 }
}
if ((opt == EUROPEANPUT ) ) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<=i;j++) {
      C[j] = dpu*C[j] + dpd*C[j+1];
      St[j]=St[j]*oneoveru;
   }
#ifdef EUROPEANOPTION_DEBUG
 PRINTMATRIX("C ","I   ",C,(i+1)); 
#endif
 }
}

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



