diff --git a/src/search.cpp b/src/search.cpp
index 7736c4121..5ec0def6a 100644
--- a/src/search.cpp
+++ b/src/search.cpp
@@ -96,6 +96,32 @@ int correction_value(const Worker& w, const Position& pos, const Stack* const ss
     return 6995 * pcv + 6593 * micv + 7753 * (wnpcv + bnpcv) + 6049 * cntcv;
 }
 
+int risk_tolerance(const Position& pos, Value v) {
+    // Returns (some constant of) second derivative of sigmoid.
+    static constexpr auto sigmoid_d2 = [](int x, int y) {
+        return -345600 * x / (x * x + 3 * y * y);
+    };
+
+    int material = pos.count<PAWN>() + 3 * pos.count<KNIGHT>() + 3 * pos.count<BISHOP>()
+                 + 5 * pos.count<ROOK>() + 9 * pos.count<QUEEN>();
+
+    int m = std::clamp(material, 17, 78);
+
+    // a and b are the crude approximation of the wdl model.
+    // The win rate is: 1/(1+exp((a-v)/b))
+    // The loss rate is 1/(1+exp((v+a)/b))
+    int a = ((-m * 3220 / 256 + 2361) * m / 256 - 586) * m / 256 + 421;
+    int b = ((m * 7761 / 256 - 2674) * m / 256 + 314) * m / 256 + 51;
+
+
+    // The risk utility is therefore d/dv^2 (1/(1+exp(-(v-a)/b)) -1/(1+exp(-(-v-a)/b)))
+    // -115200x/(x^2+3) = -345600(ab) / (a^2+3b^2) (multiplied by some constant) (second degree pade approximant)
+    int winning_risk = sigmoid_d2(v - a, b);
+    int losing_risk  = -sigmoid_d2(-v - a, b);
+
+    return (winning_risk + losing_risk) * 60 / b;
+}
+
 // Add correctionHistory value to raw staticEval and guarantee evaluation
 // does not hit the tablebase range.
 Value to_corrected_static_eval(const Value v, const int cv) {
@@ -1166,6 +1192,9 @@ moves_loop:  // When in check, search starts here
 
         r -= std::abs(correctionValue) / 31568;
 
+        if (PvNode && !is_decisive(bestValue))
+            r -= risk_tolerance(pos, bestValue);
+
         // Increase reduction for cut nodes
         if (cutNode)
             r += 2608 + 1024 * !ttData.move;