Cleanup, updated README

Makefile

@@ -3,6 +3,7 @@ CFLAGS=-Wall -g -fPIC -O3 -pthread
 CPPFLAGS=-pedantic -std=c++11
 LDFLAGS=-lm -lstdc++
 PREFIX=/usr/local
+THIS_DIR=$(dir $(abspath $(lastword $(MAKEFILE_LIST))))
 
 WINDOWS_CC=x86_64-w64-mingw32-gcc
 WINDOWS_CFLAGS= -Wall -O3 -ansi -pedantic-errors -pthread
@@ -55,8 +56,7 @@ windows: build/opt-pricer.exe
 .PHONY: install
 install : build/opt-pricer
 	@mkdir -p $(DESTDIR)$(PREFIX)/bin
-	@cp $< $(DESTDIR)$(PREFIX)/bin/opt-pricer
-	@rm -f $<
+	@ln -sf $(THIS_DIR)$< $(DESTDIR)$(PREFIX)/bin/opt-pricer
 
 .PHONY: uninstall
 uninstall :
@@ -65,5 +65,3 @@ uninstall :
 .PHONY: clean
 clean :
 	@rm -rf build
-	@mkdir build
-	@mkdir build/depends

README

@@ -2,7 +2,7 @@ Demonstration of a Black-Scholes Model in C
 
 Note that there are two models, a pure Monte Carlo implementation
 and the closed form Black-Scholes equation. Note that the code runs
-multi-threaded.
+multi-threaded and the random numbers are a sobol sequence.
 
 To install:
 	$ git clone https://github.com/kevindkeogh/opt-pricer.git
@@ -29,12 +29,12 @@ To install:
 		    |Risk-free:  |       3.00% |          3.00% |
 		    |Implied Vol:|      20.00% |         20.00% |
 		    ---------------------------------------------
-		    |Fair value: |    11.8866  |       11.8941  |
-		    |Delta:      |     0.6202  |        0.6206  |
-		    |Gamma:      |     0.0253  |        0.0134  |
-		    |Vega:       |     0.4660  |        0.4668  |
-		    |Theta:      |    -4.6138  |       -4.6235  |
-		    |Rho:        |     0.7513  |        0.7699  |
+		    |Fair value: |    11.8866  |       11.8866  |
+		    |Delta:      |     0.6202  |        0.6202  |
+		    |Gamma:      |     0.0253  |        0.0157  |
+		    |Vega:       |     0.4660  |        0.4660  |
+		    |Theta:      |    -4.6138  |       -4.6154  |
+		    |Rho:        |     0.7513  |        0.7630  |
 		    |Simulations:|             |   100,000,000  |
 		    ---------------------------------------------
 
@@ -42,7 +42,15 @@ To install:
 Recommended number of simulations is 100,000,000 for Gamma convergence, the
 other Greeks converge by 1,000,000.
 
+LICENSE:
+	My code, which comprises most of the routine functionality, is MIT.
+	There are two functions,
+		i8_sobol_generate - generate sobol sequence
+		r8_normal_01_cdf_inverse - calculate inverse normal cdf
+	which are LGPL. The Makefile compiles them as a separate shared object
+	and dynamically links them. The files are included (with minor
+	adjustments to make it compile) in src.
+
 TODO:
- 1. Find a C implementation for Sobol sequence so that Gamma convergence
-    is faster
- 2. Tests...
+	1. Tests...
+	2. Update windows cross-compilation, currently broken

src/gbm_mc.c

@@ -38,7 +38,6 @@ void *run_simulations(void *opt_ptr)
 	double gamma_shift = 0, base_gamma = 0, upper_gamma = 0, lower_gamma = 0;
 	double price, delta, gamma, vega, theta, rho;
 	double base = 0;
-	double max_rand = 0;
 
 	struct Option *opt = (struct Option *) opt_ptr;
 
@@ -52,7 +51,7 @@ void *run_simulations(void *opt_ptr)
 
 	for (i = 0; i < opt->sims; i++) {
 		rand = opt->randoms[i];
-		max_rand = rand > max_rand ? rand : max_rand;
+
 		/* Base scenario */
 		level = gbm_simulation(opt->spot, opt->rfr, opt->vol, tte, rand);
 		base += max((level - opt->strike) * opt->type, 0);
@@ -83,13 +82,13 @@ void *run_simulations(void *opt_ptr)
 		rho_shift += max((level - opt->strike) * opt->type, 0);
 	}
 
-	df = 1 / pow((1 + opt->rfr), tte);
+	df = exp(-opt->rfr * tte);
 	price = base / opt->sims * df;
 	delta = delta_shift / opt->sims * df;
 	gamma = gamma_shift / opt->sims * df;
 	vega = vega_shift / opt->sims * df;
-	theta = theta_shift / opt->sims * 1 / pow((1 + opt->rfr), tte - 1./365.);
-	rho = rho_shift / opt->sims * 1 / pow((1 + opt->rfr + 0.01), tte);
+	theta = theta_shift / opt->sims * exp(-opt->rfr * theta_tte);
+	rho = rho_shift / opt->sims * exp(-1 * (opt->rfr + 0.01)  * tte);
 
 	opt->fv = price;
 	opt->delta = (delta - price) / 0.0001;
@@ -125,39 +124,20 @@ void gbm(struct Option *opt)
 		}
 	}
 
-	/* Fill 2D array with normal randoms
-	 * The purpose is to give the Option structs the random
-	 * numbers they will need for simulations, as opposed to
-	 * having the individual threads do so (rand() is not thread-safe)
-	j = 0;
-	i = 0;
-	for(i = 0; i < NUM_THREADS;){
-		randoms[i][j] = gaussrand();
-		if (j >= (opt->sims / NUM_THREADS)) {
-			j = 0;
-			i += 1;
-		}
-		j++;
-	}
-	*/
-	/* Set the number of simulations on a per-thread basis
-	 */
+	/* Set the number of simulations on a per-thread basis */
 	opt->sims = opt->sims / NUM_THREADS;
+	threads = malloc(sizeof(pthread_t) * NUM_THREADS);
 
 	for(i=0; i<NUM_THREADS; i++) {
 		options[i] = *opt;
-		/* These are pointers, so need to copy them directly,
+		/* The dates pointers, so need to copy them directly,
 		 * otherwise we will probably have 2 threads trying to
 		 * access them at the same time
 		 */
 		options[i].expiry_date = opt->expiry_date;
 		options[i].value_date = opt->value_date;
 		options[i].randoms = randoms[i];
-	}
 
-	threads = malloc(sizeof(pthread_t) * NUM_THREADS);
-
-	for(i=0; i<NUM_THREADS; i++) {
 		if (pthread_create(&threads[i], NULL, run_simulations, &options[i])) {
 			printf("Error in thread creation\n");
 		}
@@ -180,9 +160,8 @@ void gbm(struct Option *opt)
 		opt->sims  += result->sims;
 	}
 
-	free(threads);
 	for(i=0; i<NUM_THREADS; i++)
 		free(randoms[i]);
 	free(randoms);
-
+	free(threads);
 }