diff --git a/models.go b/models.go
index 12932a5..602e837 100644
--- a/models.go
+++ b/models.go
@@ -82,7 +82,33 @@ func (opt *Option) PriceMonteCarlo() {
 	opt.Delta = (opt.Delta/float64(opt.Sims)*df - opt.FV) / 0.0001
 	opt.Gamma = (opt.Gamma/float64(opt.Sims)*df - opt.FV) / 10000
 	opt.Vega = (opt.Vega/float64(opt.Sims)*df - opt.FV) / 0.01
-	opt.Theta = (opt.Theta/float64(opt.Sims)*math.Exp(-opt.Rfr*(tte-1./365)) - opt.FV) / -(1. / 365)
+	opt.Theta = (opt.Theta/float64(opt.Sims)*math.Exp(-opt.Rfr*(tte-1./365)) -
+		opt.FV) / -(1. / 365)
 	opt.Rho = (opt.Rho/float64(opt.Sims)*math.Exp(-(opt.Rfr+0.01)*tte) - opt.FV)
 }
 
+func (opt *Option) PriceClosedForm() {
+	var d1, d2, tte float64
+
+	tte = opt.ExpiryDate.Sub(opt.ValueDate).Hours() / (24 * 365)
+
+	d1 = (math.Log(opt.Spot/opt.Strike) + tte*(opt.Rfr+math.Pow(opt.Vol, 2)/2)) /
+		(opt.Vol * math.Pow(tte, 0.5))
+	d2 = d1 - (opt.Vol * math.Pow(tte, 0.5))
+
+	opt.FV = opt.Spot*NormCDF(d1*float64(opt.OptType))*float64(opt.OptType) -
+		opt.Strike*math.Exp(-opt.Rfr*tte)*
+			NormCDF(d2*float64(opt.OptType))*float64(opt.OptType)
+	opt.Delta = NormCDF(d1)
+	if opt.OptType == -1 {
+		opt.Delta -= 1
+	}
+	opt.Gamma = 1 / (opt.Vol * math.Sqrt(tte) * math.Sqrt(2*math.Pi)) *
+		math.Exp(-d1*d1/2)
+	opt.Vega = opt.Spot / 100 * math.Pow(tte, 0.5) * NormPDF(d1)
+	opt.Rho = float64(opt.OptType) * tte * opt.Strike * math.Exp(-opt.Rfr*tte) *
+		NormCDF(float64(opt.OptType)*d2) / 100
+	opt.Theta = -float64(opt.OptType)*(opt.Rfr*opt.Strike*math.Exp(-opt.Rfr*tte)*
+		NormCDF(float64(opt.OptType)*d1)) - (opt.Vol/2*math.Pow(tte, 0.5))*
+		opt.Spot*NormPDF(float64(opt.OptType)*d1)
+}
diff --git a/serve.go b/serve.go
index d90324e..f922c27 100644
--- a/serve.go
+++ b/serve.go
@@ -22,6 +22,7 @@ func index(w http.ResponseWriter, r *http.Request) {
 		json.Unmarshal(body, &mcOpt)
 
 		mcOpt.PriceMonteCarlo()
+		bsOpt.PriceClosedForm()
 
 		options := make(map[string]Option)
 		options["MonteCarlo"] = mcOpt
diff --git a/utils.go b/utils.go
new file mode 100644
index 0000000..2f36bb1
--- /dev/null
+++ b/utils.go
@@ -0,0 +1,35 @@
+package main
+
+import "math"
+
+// NormCDF computes the Cumulative Density Function of a Standard Normal
+// Distribution.
+// Code adapted from https://www.johndcook.com/blog/cpp_phi/
+func NormCDF(z float64) float64 {
+	// constants
+	var a1, a2, a3, a4, a5, p, t, y, sign float64
+
+	a1 = 0.254829592
+	a2 = -0.284496736
+	a3 = 1.421413741
+	a4 = -1.453152027
+	a5 = 1.061405429
+	p = 0.3275911
+
+	sign = 1
+	if z < 0 {
+		sign = -1
+	}
+	x := math.Abs(z) / math.Sqrt2
+
+	t = 1.0 / (1.0 + p*x)
+	y = 1.0 - (((((a5*t+a4)*t)+a3)*t+a2)*t+a1)*t*math.Exp(-x*x)
+
+	return 0.5 * (1.0 + sign*y)
+}
+
+// NormPDF computes the Probability Density Function of a Standard Normal
+// Distribution.
+func NormPDF(z float64) float64 {
+	return 1 / math.Sqrt(2*math.Pi) * math.Exp(-math.Pow(z, 2)/2)
+}