--- a/src/actuarialmath/insurance.py
+++ b/src/actuarialmath/insurance.py
@@ -108,7 +108,7 @@
         """
         if moment == self.VARIANCE:
             A2 = self.whole_life_insurance(x, s=s, moment=2, discrete=discrete)
-            A1 = self.whole_life_insurance(x, s=s, discrete=discrete)**2
+            A1 = self.whole_life_insurance(x, s=s, discrete=discrete)
             return self.insurance_variance(A2=A2, A1=A1, b=b)
         return self.A_x(x, s=s, t=self.WHOLE, benefit=lambda x,t: b, 
                         moment=moment, discrete=discrete)
--- /dev/null
+++ b/tests/test_whole_life_variance.py
@@ -0,0 +1,20 @@
+"""Regression against an exactly enumerable two-year lifetime distribution."""
+from fractions import Fraction
+import sys
+
+import pytest
+from actuarialmath import LifeTable
+
+
+@pytest.mark.parametrize('q_text', ['0', '0.25', '0.5', '1'])
+@pytest.mark.parametrize('rate_text', ['0.01', '0.05'])
+@pytest.mark.parametrize('benefit', [1, 100])
+def test_discrete_whole_life_variance(q_text, rate_text, benefit):
+    """A payout at year one or two has the Bernoulli two-point variance."""
+    q, rate = Fraction(q_text), Fraction(rate_text)
+    v = 1 / (1 + rate)
+    expected = q * (1 - q) * (benefit * (v - v * v)) ** 2
+    life = LifeTable().set_table(q={40: float(q), 41: 1}).set_interest(i=float(rate))
+    actual = life.whole_life_insurance(40, b=benefit, moment=life.VARIANCE)
+    tolerance = 32 * sys.float_info.epsilon * benefit ** 2
+    assert actual == pytest.approx(float(expected), rel=0, abs=tolerance)