diff --git a/Makefile b/Makefile
index ce147fc..022a728 100644
--- a/Makefile
+++ b/Makefile
@@ -1,35 +1,47 @@
 CC=gcc
-WINDOWS_CC=x86_64-w64-mingw32-gcc
 CFLAGS=-Wall -fPIC -O3 -ansi -pedantic-errors
-WINDOWS_CFLAGS= -Wall -O3 -ansi -pedantic-errors
 LDFLAGS=-lm
 PREFIX= /usr/local
 
-opt-pricer : src/opt-pricer.c gbm.o black_scholes.o utils.o
+WINDOWS_CC=x86_64-w64-mingw32-gcc
+WINDOWS_CFLAGS= -Wall -O3 -ansi -pedantic-errors
+WINDOWS_PLATFORM= windows
+
+.DEFAULT_GOAL := build/opt-pricer
+
+build/opt-pricer : src/opt-pricer.c build/depends/linux/gbm.o build/depends/linux/black_scholes.o build/depends/linux/utils.o
 	@$(CC) $(CFLAGS) $^ $(LDFLAGS) -o $@
 
-gbm.o : src/gbm_mc.c
+build/depends/%/gbm.o : src/gbm_mc.c | folders
 	@$(CC) $(CFLAGS) -c $^ $(LDFLAGS) -o $@
 
-black_scholes.o : src/black_scholes.c
+build/depends/%/black_scholes.o : src/black_scholes.c | folders
 	@$(CC) $(CFLAGS) -c $^ $(LDFLAGS) -o $@
 
-utils.o : src/utils.c
+build/depends/%/utils.o : src/utils.c | folders
 	@$(CC) $(CFLAGS) -c $^ $(LDFLAGS) -o $@
 
-strptime.o : src/strptime.c
+build/depends/%/strptime.o : src/strptime.c | folders
 	@$(CC) $(CFLAGS) -c $^ $(LDFLAGS) -o $@
 
-opt-pricer.exe : CC=$(WINDOWS_CC)
-opt-pricer.exe : CFLAGS=$(WINDOWS_CFLAGS)
-opt-pricer.exe : src/opt-pricer.c gbm.o black_scholes.o utils.o strptime.o
-	@$(WINDOWS_CC) $(WINDOWS_CFLAGS) $^ $(LDFLAGS) -o $@
+folders:
+	@mkdir -p build
+	@mkdir -p build/depends
+	@mkdir -p build/depends/linux
+	@mkdir -p build/depends/windows
+
+
+build/opt-pricer.exe : CC=$(WINDOWS_CC)
+build/opt-pricer.exe : CFLAGS=$(WINDOWS_CFLAGS)
+build/opt-pricer.exe : src/opt-pricer.c build/depends/windows/gbm.o build/depends/windows/black_scholes.o build/depends/windows/utils.o build/depends/windows/strptime.o
+	@$(CC) $(CFLAGS) $^ $(LDFLAGS) -o $@
+
 
 .PHONY: windows
-windows: opt-pricer.exe
+windows: build/opt-pricer.exe
 
 .PHONY: install
-install : opt-pricer
+install : build/opt-pricer
 	@mkdir -p $(DESTDIR)$(PREFIX)/bin
 	@cp $< $(DESTDIR)$(PREFIX)/bin/opt-pricer
 	@rm -f $<
@@ -40,4 +52,6 @@ uninstall :
 
 .PHONY: clean
 clean :
-	@rm -f gbm.o black_scholes.o utils.o strptime.o
+	@rm -rf build
+	@mkdir build
+	@mkdir build/depends
diff --git a/README b/README
new file mode 100644
index 0000000..2e3b528
--- /dev/null
+++ b/README
@@ -0,0 +1,12 @@
+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.
+
+Recommended number of simulations is 100,000,000 for Gamma convergence, the
+other Greeks converge by 1,000,000.
+
+TODO:
+ 1. Find a C implementation for Sobol sequence so that Gamma convergence
+    is faster
+ 2. Tests...
diff --git a/src/black_scholes.c b/src/black_scholes.c
index 6b3c9e5..cd56c40 100644
--- a/src/black_scholes.c
+++ b/src/black_scholes.c
@@ -22,5 +22,9 @@ void bsm(struct Option *opt)
 	price = opt->spot * normalcdf(d1 * opt->type) * opt->type - opt->strike *
 		exp(-opt->rfr * tte) * normalcdf(d2 * opt->type) * opt->type;
 	opt->fv = price;
-	opt->vega = opt->spot / 100 * exp(-opt->rfr * tte) * pow(tte, 0.5) * normalpdf(d1);
+	opt->delta = (opt->type == 1) ? normalcdf(d1) : normalcdf(d1) - 1;
+	opt->gamma = 1 / (opt->spot * opt->vol * pow(tte, 0.5)) * normalcdf(d1);
+	opt->vega = opt->spot / 100 * pow(tte, 0.5) * normalpdf(d1);
+	opt->rho = opt->type * tte * opt->strike * exp(-opt->rfr * tte) * normalcdf(opt->type * d2) / 100;
+	opt->theta = -opt->type * (opt->rfr * opt->strike * exp(-opt->rfr * tte) * normalcdf(opt->type * d2)) - (opt->vol / (2 * pow(tte, 0.5)) * opt->spot * normalpdf(opt->type * d1));
 }
diff --git a/src/gbm_mc.c b/src/gbm_mc.c
index 46bcd72..b2224a0 100644
--- a/src/gbm_mc.c
+++ b/src/gbm_mc.c
@@ -19,9 +19,10 @@ double gbm_simulation(double spot, double rfr, double vol, double tte, double ra
 
 void gbm(struct Option *opt)
 {
-	double tte, expiry_date, value_date, level, rand;
+	double tte, theta_tte, expiry_date, value_date, level, rand, df;
 	double delta_shift = 0, vega_shift = 0, theta_shift = 0, rho_shift = 0;
-	double price, delta, vega, theta, rho;
+	double gamma_shift = 0, base_gamma = 0, upper_gamma = 0, lower_gamma = 0;
+	double price, delta, gamma, vega, theta, rho;
 	double base = 0;
 	int i;
 
@@ -30,30 +31,52 @@ void gbm(struct Option *opt)
 	expiry_date = mktime(opt->expiry_date);
 	value_date = mktime(opt->value_date);
 	tte = difftime(expiry_date, value_date) / (60 * 60 * 24 * 365);
+	theta_tte = tte - 1. / 365;
 
 	for (i=0; i<opt->sims; i++) {
 		rand = gaussrand();
+		/* Base scenario */
 		level = gbm_simulation(opt->spot, opt->rfr, opt->vol, tte, rand);
 		base += max((level - opt->strike) * opt->type, 0);
-		level = gbm_simulation(opt->spot + 1, opt->rfr, opt->vol, tte, rand);
+
+		/* Delta scenario */
+		level = gbm_simulation(opt->spot + 0.0001, opt->rfr, opt->vol, tte, rand);
 		delta_shift += max((level - opt->strike) * opt->type, 0);
-		level = gbm_simulation(opt->spot, opt->rfr, opt->vol + 0.01, tte, rand);
+
+		/* Gamma scenario */
+		upper_gamma = gbm_simulation(opt->spot + 0.00001, opt->rfr, opt->vol, tte, rand);
+		base_gamma = gbm_simulation(opt->spot, opt->rfr, opt->vol, tte, rand);
+		lower_gamma = gbm_simulation(opt->spot - 0.00001, opt->rfr, opt->vol, tte, rand);
+		gamma_shift += (max((upper_gamma - opt->strike) * opt->type, 0) -
+			2 * max((base_gamma - opt->strike) * opt->type, 0) +
+			max((lower_gamma - opt->strike) * opt->type, 0)) /
+			pow(0.00001, 2);
+
+		/* Vega scenario */
+		level = gbm_simulation(opt->spot, opt->rfr, opt->vol + 0.0001, tte, rand);
 		vega_shift += max((level - opt->strike) * opt->type, 0);
-		level = gbm_simulation(opt->spot, opt->rfr, opt->vol, tte - 1./365., rand);
+
+		/* Theta scenario */
+		level = gbm_simulation(opt->spot, opt->rfr, opt->vol, theta_tte, rand);
 		theta_shift += max((level - opt->strike) * opt->type, 0);
+
+		/* Rho scenario */
 		level = gbm_simulation(opt->spot, opt->rfr + 0.01, opt->vol, tte, rand);
 		rho_shift += max((level - opt->strike) * opt->type, 0);
 	}
 
-	price = base / opt->sims * 1 / pow((1 + opt->rfr), tte);
-	delta = delta_shift / opt->sims * 1 / pow((1 + opt->rfr), tte);
-	vega = vega_shift / opt->sims * 1 / pow((1 + opt->rfr), tte);
+	df = 1 / pow((1 + 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);
 
 	opt->fv = price;
-	opt->delta = (delta - price);
-	opt->vega = (vega - price);
-	opt->theta = (theta - price);
+	opt->delta = (delta - price) / 0.0001;
+	opt->gamma = gamma;
+	opt->vega = (vega - price) / 0.01;
+	opt->theta = (theta - price) / (tte - theta_tte);
 	opt->rho = (rho - price);
 }
diff --git a/src/opt-pricer.c b/src/opt-pricer.c
index 6297ffd..d3ef5f5 100644
--- a/src/opt-pricer.c
+++ b/src/opt-pricer.c
@@ -8,6 +8,7 @@
 #include "strptime.h"
 #endif
 
+#include <locale.h>
 #include <stdio.h>
 #include <string.h>
 #include <time.h>
@@ -70,6 +71,12 @@ int print_help(void)
 }
 
 
+char* read_type(struct Option opt)
+{
+	return (opt.type == 1) ? "Call" : "Put";
+}
+
+
 int main(int argc, char *argv[])
 {
 	double spot = 0, strike = 0, rfr = 0, vol = 0, sims = 1000;
@@ -158,41 +165,42 @@ int main(int argc, char *argv[])
 	strftime(value_date, 11, "%Y-%m-%d", &value);
 	bsm(&bs_opt);
 	gbm(&mc_opt);
+
+	setlocale(LC_ALL,"");
 	printf(
 		"\nValuation date: %s\n\n"
 		"	              | BS Analytic | BS Monte Carlo |\n"
 		"	 ---------------------------------------------\n"
-		"	 |Type:       | %10.1i  | %13.1i  |\n"
+		"	 |Type:       | %10s  | %13s  |\n"
 		"	 |Spot:       | %10.2f  | %13.2f  |\n"
+		"	 |Expiry:     | %s  |    %s  |\n"
 		"	 |Strike:     | %10.2f  | %13.2f  |\n"
 		"	 |Risk-free:  | %10.2f%% | %13.2f%% |\n"
 		"	 |Implied Vol:| %10.2f%% | %13.2f%% |\n"
-		/*
-		"	 |sims:     | %4.2f| %4.2f|\n"
-		"	 |exp:      | %4.2s| %4.2f|\n"
-		*/
+		"	 ---------------------------------------------\n"
 		"	 |Fair value: |   %8.4f  |   %11.4f  |\n"
 		"	 |Delta:      |   %8.4f  |   %11.4f  |\n"
-		"	 |Vega:       |   %8.4f  |   %11.4f  |\n"
-		"	 |Theta:      |   %8.4f  |   %11.4f  |\n"
-		"	 |Rho:        |   %8.4f  |   %11.4f  |\n",
+		"	 |Gamma:      |   %8.4f  |   %11.4f  |\n",
 		value_date,
-		bs_opt.type, mc_opt.type,
+		read_type(bs_opt), read_type(mc_opt),
 		bs_opt.spot, mc_opt.spot,
+		expiry_date, expiry_date,
 		bs_opt.strike, mc_opt.strike,
 		bs_opt.rfr*100, mc_opt.rfr*100,
 		bs_opt.vol*100, mc_opt.vol*100,
 		bs_opt.fv, mc_opt.fv,
 		bs_opt.delta, mc_opt.delta,
+		bs_opt.gamma, mc_opt.gamma);
+	printf(
+		"	 |Vega:       |   %8.4f  |   %11.4f  |\n"
+		"	 |Theta:      |   %8.4f  |   %11.4f  |\n"
+		"	 |Rho:        |   %8.4f  |   %11.4f  |\n"
+		"	 |Simulations:|             |   %'11ld  |\n"
+		"	 ---------------------------------------------\n\n",
 		bs_opt.vega, mc_opt.vega,
 		bs_opt.theta, mc_opt.theta,
-		bs_opt.rho, mc_opt.rho);
-	printf(
-		"	 ---------------------------------------------\n\n"
-	      );
-	/*
-	printf("%s", buffer1);
-	printf("%s", buffer2);
-	*/
+		bs_opt.rho, mc_opt.rho,
+		            mc_opt.sims);
+
 	return 0;
 }