# include # include # include # include # include "black_scholes.h" /******************************************************************************/ double *asset_path ( double s0, double mu, double sigma, double t1, int n, int *seed ) /******************************************************************************/ /* Purpose: ASSET_PATH simulates the behavior of an asset price over time. Licensing: This code is distributed under the MIT license. Modified: 18 February 2012 Author: Original MATLAB version by Desmond Higham. C version by John Burkardt. Reference: Desmond Higham, Black-Scholes for Scientific Computing Students, Computing in Science and Engineering, November/December 2004, Volume 6, Number 6, pages 72-79. Parameters: Input, double S0, the asset price at time 0. Input, double MU, the expected growth rate. Input, double SIGMA, the volatility of the asset. Input, double T1, the expiry date. Input, integer N, the number of steps to take between 0 and T1. Input/output, int *SEED, a seed for the random number generator. Output, double ASSET_PATH[N+1], the option values from time 0 to T1 in equal steps. */ { double dt; int i; double p; double *r; double *s; dt = t1 / ( double ) ( n ); r = r8vec_normal_01_new ( n, seed ); s = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) ); s[0] = s0; p = s0; for ( i = 1; i <= n; i++ ) { p = p * exp ( ( mu - sigma * sigma ) * dt + sigma * sqrt ( dt ) * r[i-1] ); s[i] = p; } free ( r ); return s; } /******************************************************************************/ double binomial ( double s0, double e, double r, double sigma, double t1, int m ) /******************************************************************************/ /* Purpose: BINOMIAL uses the binomial method for a European call. Licensing: This code is distributed under the MIT license. Modified: 18 February 2012 Author: Original MATLAB version by Desmond Higham. C version by John Burkardt. Reference: Desmond Higham, Black-Scholes for Scientific Computing Students, Computing in Science and Engineering, November/December 2004, Volume 6, Number 6, pages 72-79. Parameters: Input, double S0, the asset price at time 0. Input, double E, the exercise price. Input, double R, the interest rate. Input, double SIGMA, the volatility of the asset. Input, double T1, the expiry date. Input, int M, the number of steps to take between 0 and T1. Output, double BIONOMIAL, the option value at time 0. */ { double a; double c; double d; double dt; int i; int n; double p; double u; double *w; /* Time stepsize. */ dt = t1 / ( double ) ( m ); a = 0.5 * ( exp ( - r * dt ) + exp ( ( r + sigma * sigma ) * dt ) ); d = a - sqrt ( a * a - 1.0 ); u = a + sqrt ( a * a - 1.0 ); p = ( exp ( r * dt ) - d ) / ( u - d ); w = ( double * ) malloc ( ( m + 1 ) * sizeof ( double ) ); for ( i = 0; i <= m; i++ ) { w[i] = fmax ( s0 * pow ( d, m - i ) * pow ( u, i ) - e, 0.0 ); } /* Trace backwards to get the option value at time 0. */ for ( n = m - 1; 0 <= n; n-- ) { for ( i = 0; i <= n; i++ ) { w[i] = ( 1.0 - p ) * w[i] + p * w[i+1]; } } for ( i = 0; i < m + 1; i++ ) { w[i] = exp ( - r * t1 ) * w[i]; } c = w[0]; free ( w ); return c; } /******************************************************************************/ double bsf ( double s0, double t0, double e, double r, double sigma, double t1 ) /******************************************************************************/ /* Purpose: BSF evaluates the Black-Scholes formula for a European call. Licensing: This code is distributed under the MIT license. Modified: 18 February 2012 Author: Original MATLAB version by Desmond Higham. C version by John Burkardt. Reference: Desmond Higham, Black-Scholes for Scientific Computing Students, Computing in Science and Engineering, November/December 2004, Volume 6, Number 6, pages 72-79. Parameters: Input, double S0, the asset price at time T0. Input, double T0, the time at which the asset price is known. Input, double E, the exercise price. Input, double R, the interest rate. Input, double SIGMA, the volatility of the asset. Input, double T1, the expiry date. Output, double BSF, the value of the call option. */ { double c; double d1; double d2; double n1; double n2; double tau; tau = t1 - t0; if ( 0.0 < tau ) { d1 = ( log ( s0 / e ) + ( r + 0.5 * sigma * sigma ) * tau ) / ( sigma * sqrt ( tau ) ); d2 = d1 - sigma * sqrt ( tau ); n1 = 0.5 * ( 1.0 + erf ( d1 / sqrt ( 2.0 ) ) ); n2 = 0.5 * ( 1.0 + erf ( d2 / sqrt ( 2.0 ) ) ); c = s0 * n1 - e * exp ( - r * tau ) * n2; } else { c = fmax ( s0 - e, 0.0 ); } return c; } /******************************************************************************/ double *forward ( double e, double r, double sigma, double t1, int nx, int nt, double smax ) /******************************************************************************/ /* Purpose: FORWARD uses the forward difference method to value a European call option. Licensing: This code is distributed under the MIT license. Modified: 18 February 2012 Author: Original MATLAB version by Desmond Higham. C version by John Burkardt. Reference: Desmond Higham, Black-Scholes for Scientific Computing Students, Computing in Science and Engineering, November/December 2004, Volume 6, Number 6, pages 72-79. Parameters: Input, double E, the exercise price. Input, double R, the interest rate. Input, double SIGMA, the volatility of the asset. Input, double T1, the expiry date. Input, int NX, the number of "space" steps used to divide the interval [0,L]. Input, int NT, the number of time steps. Input, double SMAX, the maximum value of S to consider. Output, double U[(NX-1)*(NT+1)], the value of the European call option. */ { double *a; double *b; double *c; double dt; double dx; int i; int j; double p; double t; double *u; double u0; dt = t1 / ( double ) ( nt ); dx = smax / ( double ) ( nx ); a = ( double * ) malloc ( ( nx - 1 ) * sizeof ( double ) ); b = ( double * ) malloc ( ( nx - 1 ) * sizeof ( double ) ); c = ( double * ) malloc ( ( nx - 1 ) * sizeof ( double ) ); for ( i = 0; i < nx - 1; i++ ) { b[i] = 1.0 - r * dt - dt * pow ( sigma * ( i + 1 ), 2 ); } for ( i = 0; i < nx - 2; i++ ) { c[i] = 0.5 * dt * pow ( sigma * ( i + 1 ), 2 ) + 0.5 * dt * r * ( i + 1 ); } for ( i = 1; i < nx - 1; i++ ) { a[i] = 0.5 * dt * pow ( sigma * ( i + 1 ), 2 ) - 0.5 * dt * r * ( i + 1 ); } u = ( double * ) malloc ( ( nx - 1 ) * ( nt + 1 ) * sizeof ( double ) ); u0 = 0.0; for ( i = 0; i < nx - 1; i++ ) { u0 = u0 + dx; u[i+0*(nx-1)] = fmax ( u0 - e, 0.0 ); } for ( j = 0; j < nt; j++ ) { t = ( double ) ( j ) * t1 / ( double ) ( nt ); p = 0.5 * dt * ( nx - 1 ) * ( sigma * sigma * ( nx - 1 ) + r ) * ( smax - e * exp ( - r * t ) ); for ( i = 0; i < nx - 1; i++ ) { u[i+(j+1)*(nx-1)] = b[i] * u[i+j*(nx-1)]; } for ( i = 0; i < nx - 2; i++ ) { u[i+(j+1)*(nx-1)] = u[i+(j+1)*(nx-1)] + c[i] * u[i+1+j*(nx-1)]; } for ( i = 1; i < nx - 1; i++ ) { u[i+(j+1)*(nx-1)] = u[i+(j+1)*(nx-1)] + a[i] * u[i-1+j*(nx-1)]; } u[nx-2+(j+1)*(nx-1)] = u[nx-2+(j+1)*(nx-1)] + p; } free ( a ); free ( b ); free ( c ); return u; } /******************************************************************************/ double *mc ( double s0, double e, double r, double sigma, double t1, int m, int *seed ) /******************************************************************************/ /* Purpose: MC uses Monte Carlo valuation on a European call. Licensing: This code is distributed under the MIT license. Modified: 18 February 2012 Author: Original MATLAB version by Desmond Higham. C version by John Burkardt. Reference: Desmond Higham, Black-Scholes for Scientific Computing Students, Computing in Science and Engineering, November/December 2004, Volume 6, Number 6, pages 72-79. Parameters: Input, double S0, the asset price at time 0. Input, double E, the exercise price. Input, double R, the interest rate. Input, double SIGMA, the volatility of the asset. Input, double T1, the expiry date. Input, int M, the number of simulations. Input/output, int *SEED, a seed for the random number generator. Output, double MC[2], the estimated range of the valuation. */ { double *conf; int i; double pmean; double *pvals; double std; double *svals; double *u; double width; u = r8vec_normal_01_new ( m, seed ); svals = ( double * ) malloc ( m * sizeof ( double ) ); for ( i = 0; i < m; i++ ) { svals[i] = s0 * exp ( ( r - 0.5 * sigma * sigma ) * t1 + sigma * sqrt ( t1 ) * u[i] ); } pvals = ( double * ) malloc ( m * sizeof ( double ) ); for ( i = 0; i < m; i++ ) { pvals[i] = exp ( - r * t1 ) * fmax ( svals[i] - e, 0.0 ); } pmean = 0.0; for ( i = 0; i < m; i++ ) { pmean = pmean + pvals[i]; } pmean = pmean / ( double ) ( m ); std = 0.0; for ( i = 0; i < m; i++ ) { std = std + pow ( pvals[i] - pmean, 2 ); } std = sqrt ( std / ( double ) ( m - 1 ) ); width = 1.96 * std / sqrt ( ( double ) ( m ) ); conf = ( double * ) malloc ( 2 * sizeof ( double ) ); conf[0] = pmean - width; conf[1] = pmean + width; free ( pvals ); free ( svals ); free ( u ); return conf; } /******************************************************************************/ double *r8vec_normal_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: R8VEC_NORMAL_01_NEW returns a unit pseudonormal R8VEC. Discussion: The standard normal probability distribution function (PDF) has mean 0 and standard deviation 1. This routine can generate a vector of values on one call. It has the feature that it should provide the same results in the same order no matter how we break up the task. Before calling this routine, the user may call RANDOM_SEED in order to set the seed of the random number generator. The Box-Muller method is used, which is efficient, but generates an even number of values each time. On any call to this routine, an even number of new values are generated. Depending on the situation, one value may be left over. In that case, it is saved for the next call. Licensing: This code is distributed under the MIT license. Modified: 18 February 2012 Author: John Burkardt Parameters: Input, int N, the number of values desired. If N is negative, then the code will flush its internal memory; in particular, if there is a saved value to be used on the next call, it is instead discarded. This is useful if the user has reset the random number seed, for instance. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_NORMAL_01_NEW[N], a sample of the standard normal PDF. Local parameters: Local, int MADE, records the number of values that have been computed. On input with negative N, this value overwrites the return value of N, so the user can get an accounting of how much work has been done. Local, double R[N+1], is used to store some uniform random values. Its dimension is N+1, but really it is only needed to be the smallest even number greater than or equal to N. Local, int SAVED, is 0 or 1 depending on whether there is a single saved value left over from the previous call. Local, int X_LO, X_HI, records the range of entries of X that we need to compute. This starts off as 1:N, but is adjusted if we have a saved value that can be immediately stored in X(1), and so on. Local, double Y, the value saved from the previous call, if SAVED is 1. */ { # define R8_PI 3.141592653589793 int i; int m; static int made = 0; double *r; static int saved = 0; double *x; int x_hi; int x_lo; static double y = 0.0; x = ( double * ) malloc ( n * sizeof ( double ) ); /* I'd like to allow the user to reset the internal data. But this won't work properly if we have a saved value Y. I'm making a crock option that allows the user to signal explicitly that any internal memory should be flushed, by passing in a negative value for N. */ if ( n < 0 ) { made = 0; saved = 0; y = 0.0; return NULL; } else if ( n == 0 ) { return NULL; } /* Record the range of X we need to fill in. */ x_lo = 1; x_hi = n; /* Use up the old value, if we have it. */ if ( saved == 1 ) { x[0] = y; saved = 0; x_lo = 2; } /* Maybe we don't need any more values. */ if ( x_hi - x_lo + 1 == 0 ) { } /* If we need just one new value, do that here to avoid null arrays. */ else if ( x_hi - x_lo + 1 == 1 ) { r = r8vec_uniform_01_new ( 2, seed ); x[x_hi-1] = sqrt ( - 2.0 * log ( r[0] ) ) * cos ( 2.0 * R8_PI * r[1] ); y = sqrt ( - 2.0 * log ( r[0] ) ) * sin ( 2.0 * R8_PI * r[1] ); saved = 1; made = made + 2; free ( r ); } /* If we require an even number of values, that's easy. */ else if ( ( x_hi - x_lo + 1 ) % 2 == 0 ) { m = ( x_hi - x_lo + 1 ) / 2; r = r8vec_uniform_01_new ( 2*m, seed ); for ( i = 0; i <= 2*m-2; i = i + 2 ) { x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * R8_PI * r[i+1] ); x[x_lo+i ] = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * R8_PI * r[i+1] ); } made = made + x_hi - x_lo + 1; free ( r ); } /* If we require an odd number of values, we generate an even number, and handle the last pair specially, storing one in X(N), and saving the other for later. */ else { x_hi = x_hi - 1; m = ( x_hi - x_lo + 1 ) / 2 + 1; r = r8vec_uniform_01_new ( 2*m, seed ); for ( i = 0; i <= 2*m-4; i = i + 2 ) { x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * R8_PI * r[i+1] ); x[x_lo+i ] = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * R8_PI * r[i+1] ); } i = 2*m - 2; x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * R8_PI * r[i+1] ); y = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * R8_PI * r[i+1] ); saved = 1; made = made + x_hi - x_lo + 2; free ( r ); } return x; # undef R8_PI } /******************************************************************************/ void r8vec_print_part ( int n, double a[], int max_print, char *title ) /******************************************************************************/ /* Purpose: R8VEC_PRINT_PART prints "part" of an R8VEC. Discussion: The user specifies MAX_PRINT, the maximum number of lines to print. If N, the size of the vector, is no more than MAX_PRINT, then the entire vector is printed, one entry per line. Otherwise, if possible, the first MAX_PRINT-2 entries are printed, followed by a line of periods suggesting an omission, and the last entry. Licensing: This code is distributed under the MIT license. Modified: 27 February 2010 Author: John Burkardt Parameters: Input, int N, the number of entries of the vector. Input, double A[N], the vector to be printed. Input, int MAX_PRINT, the maximum number of lines to print. Input, char *TITLE, a title. */ { int i; if ( max_print <= 0 ) { return; } if ( n <= 0 ) { return; } fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); if ( n <= max_print ) { for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %8d: %14f\n", i, a[i] ); } } else if ( 3 <= max_print ) { for ( i = 0; i < max_print - 2; i++ ) { fprintf ( stdout, " %8d: %14f\n", i, a[i] ); } fprintf ( stdout, " ...... ..............\n" ); i = n - 1; fprintf ( stdout, " %8d: %14f\n", i, a[i] ); } else { for ( i = 0; i < max_print - 1; i++ ) { fprintf ( stdout, " %8d: %14f\n", i, a[i] ); } i = max_print - 1; fprintf ( stdout, " %8d: %14f ...more entries...\n", i, a[i] ); } return; } /******************************************************************************/ double *r8vec_uniform_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: R8VEC_UNIFORM_01_NEW returns a unit pseudorandom R8VEC. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the MIT license. Modified: 19 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Second Edition, Springer, 1987, ISBN: 0387964673, LC: QA76.9.C65.B73. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, December 1986, pages 362-376. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation, edited by Jerry Banks, Wiley, 1998, ISBN: 0471134031, LC: T57.62.H37. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, Number 2, 1969, pages 136-143. Parameters: Input, int N, the number of entries in the vector. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values. */ { int i; const int i4_huge = 2147483647; int k; double *r; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8VEC_UNIFORM_01_NEW - Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r[i] = ( double ) ( *seed ) * 4.656612875E-10; } return r; } /******************************************************************************/ void r8vec_write ( char *output_filename, int n, double x[] ) /******************************************************************************/ /* Purpose: R8VEC_WRITE writes an R8VEC file. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 10 July 2011 Author: John Burkardt Parameters: Input, char *OUTPUT_FILENAME, the output filename. Input, int N, the number of points. Input, double X[N], the data. */ { int j; FILE *output; /* Open the file. */ output = fopen ( output_filename, "wt" ); if ( !output ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8VEC_WRITE - Fatal error!\n" ); fprintf ( stderr, " Could not open the output file.\n" ); exit ( 1 ); } /* Write the data. */ for ( j = 0; j < n; j++ ) { fprintf ( output, " %24.16g\n", x[j] ); } /* Close the file. */ fclose ( output ); return; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }