diff --git a/financepy/models/black_scholes_analytic.py b/financepy/models/black_scholes_analytic.py
index 8632626e..253713a3 100644
--- a/financepy/models/black_scholes_analytic.py
+++ b/financepy/models/black_scholes_analytic.py
@@ -284,7 +284,7 @@ def bs_vanna(
     kk = k * np.exp(-r * t)
     d1 = np.log(ss / kk) / v_sqrt_t + v_sqrt_t / 2.0
     d2 = d1 - v_sqrt_t
-    vanna = np.exp(-q * t) * sqrt_t * normcdf_prime_vect(d1) * (d2 / v)
+    vanna = -np.exp(-q * t) * normcdf_prime_vect(d1) * (d2 / v)
     return vanna
 
 
diff --git a/regression_tests/golden/TestFinEquityVanillaOption_GOLDEN.testLog b/regression_tests/golden/TestFinEquityVanillaOption_GOLDEN.testLog
index f76b5282..8e6c3720 100644
--- a/regression_tests/golden/TestFinEquityVanillaOption_GOLDEN.testLog
+++ b/regression_tests/golden/TestFinEquityVanillaOption_GOLDEN.testLog
@@ -14,15 +14,15 @@ RESULTS,100000,12.50152120,12.48345880,12.48846727,0.00651217,
 RESULTS,100000,7.34775277,7.32039491,7.33606550,0.00662732,
 RESULTS,100000,4.01048407,3.98618166,3.99511073,0.00656462,
 HEADER,STOCK PRICE,CALL_VALUE_BS,CALL_DELTA_BS,CALL_VEGA_BS,CALL_THETA_BS,CALL_RHO_BS,CALL_VANNA_BS,
-RESULTS,80,1.65961061,0.19482665,15.49364885,-5.22707951,6.90602823,-0.68946422,
-RESULTS,90,4.50529071,0.38049989,24.05757490,-8.42161920,14.74763192,-0.45483484,
-RESULTS,100,9.30205599,0.57620826,27.40342752,-10.12888855,23.96081450,-0.01072075,
-RESULTS,110,15.91532099,0.73870008,24.88320231,-9.98133932,32.40231651,0.33133703,
+RESULTS,80,1.65961061,0.19482665,15.49364885,-5.22707951,6.90602823,0.97908158,
+RESULTS,90,4.50529071,0.38049989,24.05757490,-8.42161920,14.74763192,0.64589343,
+RESULTS,100,9.30205599,0.57620826,27.40342752,-10.12888855,23.96081450,0.01522413,
+RESULTS,110,15.91532099,0.73870008,24.88320231,-9.98133932,32.40231651,-0.47051895,
 HEADER,STOCK PRICE,PUT_VALUE_BS,PUT_DELTA_BS,PUT_VEGA_BS,PUT_THETA_BS,PUT_RHO_BS,PUT_VANNA_BS,
-RESULTS,80,19.60637481,-0.80022673,15.49364885,-1.14557052,-41.46859405,-0.68946422,
-RESULTS,90,12.50152120,-0.61455348,24.05757490,-4.43961554,-33.62699036,-0.45483484,
-RESULTS,100,7.34775277,-0.41884512,27.40342752,-6.24639022,-24.41380778,-0.01072075,
-RESULTS,110,4.01048407,-0.25635329,24.88320231,-6.19834633,-15.97230577,0.33133703,
+RESULTS,80,19.60637481,-0.80022673,15.49364885,-1.14557052,-41.46859405,0.97908158,
+RESULTS,90,12.50152120,-0.61455348,24.05757490,-4.43961554,-33.62699036,0.64589343,
+RESULTS,100,7.34775277,-0.41884512,27.40342752,-6.24639022,-24.41380778,0.01522413,
+RESULTS,110,4.01048407,-0.25635329,24.88320231,-6.19834633,-15.97230577,-0.47051895,
 HEADER,OPT_TYPE,TEXP,STOCK_PRICE,STRIKE,INTRINSIC,VALUE,INPUT_VOL,IMPLIED_VOL,
 RESULTS,OptionTypes.EUROPEAN_CALL,0.00300000,100.00000000,50.00000000,49.99863001,49.99863001,0.01000000,-999,
 RESULTS,OptionTypes.EUROPEAN_CALL,0.00300000,100.00000000,50.00000000,100.00000000,49.99863001,0.01000000,-999,
diff --git a/unit_tests/test_FinEquityVanillaOption.py b/unit_tests/test_FinEquityVanillaOption.py
index 48e515fd..2f7b4f8f 100644
--- a/unit_tests/test_FinEquityVanillaOption.py
+++ b/unit_tests/test_FinEquityVanillaOption.py
@@ -1,5 +1,7 @@
 # Copyright (C) 2018, 2019, 2020 Dominic O'Kane
 
+import math
+
 from financepy.utils.global_types import OptionTypes
 from financepy.products.equity.equity_vanilla_option import EquityVanillaOption
 from financepy.market.curves.discount_curve_flat import DiscountCurveFlat
@@ -73,3 +75,47 @@ def test_put_option():
 
     v = put_option.value(value_date, stock_price, discount_curve, dividend_curve, model)
     assert v.round(4) == 7.3478
+
+
+########################################################################################
+
+
+def closed_form_vanna(s, k, t, r, q, v):
+    # Black-Scholes vanna = d(vega)/d(spot) = -exp(-q t) n(d1) d2 / v
+    d1 = (math.log(s / k) + (r - q + 0.5 * v * v) * t) / (v * math.sqrt(t))
+    d2 = d1 - v * math.sqrt(t)
+    n_d1 = math.exp(-0.5 * d1 * d1) / math.sqrt(2.0 * math.pi)
+    return -math.exp(-q * t) * n_d1 * d2 / v
+
+
+def fd_vanna(option, s):
+    # d(delta)/d(vol) by central difference of the library's own delta
+    h = 1e-4
+    up = option.delta(
+        value_date, s, discount_curve, dividend_curve, BlackScholes(volatility + h)
+    )
+    dn = option.delta(
+        value_date, s, discount_curve, dividend_curve, BlackScholes(volatility - h)
+    )
+    return (up - dn) / (2 * h)
+
+
+def test_vanna():
+
+    t = (expiry_date - value_date) / 365.0
+
+    # OTM / ATM / ITM: vanna is positive below the strike and negative above it
+    for s, expected in ((80, 0.9791), (100, 0.0152), (110, -0.4705)):
+
+        call = EquityVanillaOption(expiry_date, 100.0, OptionTypes.EUROPEAN_CALL)
+        put = EquityVanillaOption(expiry_date, 100.0, OptionTypes.EUROPEAN_PUT)
+
+        v_call = float(call.vanna(value_date, s, discount_curve, dividend_curve, model))
+        v_put = float(put.vanna(value_date, s, discount_curve, dividend_curve, model))
+
+        ref = closed_form_vanna(s, 100.0, t, interest_rate, dividend_yield, volatility)
+
+        assert round(v_call, 4) == expected
+        assert abs(v_call - ref) < 1e-6
+        assert abs(v_call - fd_vanna(call, s)) < 1e-3
+        assert abs(v_call - v_put) < 1e-9  # vanna does not depend on option type