QuantLib-Python: Hull-White one-factor model calibration

This Python program is presenting the process of calibrating Hull-White One-factor interest rate model to a given set of Swaption volatilities. In the example program, I have used exactly the same data as used in the book QuantLib Python Cookbook by G. Balaraman and L. Ballabio in the Chapter 14 (Short Interest Rate Model Calibration) and corresponding results are perfectly replicated here.

Thanks for reading my blog.
-Mike

from QuantLib import *
import numpy as Numpy

# hard-coded data
def CreateSwaptionVolatilityList():
    vol = []
    vol.append(0.1148)
    vol.append(0.1108)
    vol.append(0.1070)
    vol.append(0.1021)
    vol.append(0.1000)
    return vol

class ModelCalibrator:
    def __init__(self, endCriteria):        
        self.endCriteria = endCriteria
        self.helpers = []

    def AddCalibrationHelper(self, helper):
        self.helpers.append(helper)

    def Calibrate(self, model, engine, curve, fixedParameters):
        # assign pricing engine to all calibration helpers
        for i in range(len(self.helpers)):
            self.helpers[i].setPricingEngine(engine)
        method = LevenbergMarquardt()
        if(len(fixedParameters) == 0):
            model.calibrate(self.helpers, method, self.endCriteria)
        else:
            model.calibrate(self.helpers, method, self.endCriteria,
                NoConstraint(), [], fixedParameters)

# general parameters
tradeDate = Date(15, February, 2002)
Settings.instance().evaluationDate = tradeDate
settlementDate = Date(19, February, 2002)
calendar = TARGET()
dayCounter = Actual360()

# create market data: term structure and diagonal volatilities
curveHandle = YieldTermStructureHandle(FlatForward(settlementDate, 0.04875825, Actual365Fixed()))
vol = CreateSwaptionVolatilityList()

# create calibrator object
endCriteria = EndCriteria(10000, 100, 0.000001, 0.00000001, 0.00000001)
calibrator = ModelCalibrator(endCriteria)

# add swaption helpers to calibrator
for i in range(len(vol)):
    t = i + 1; tenor = len(vol) - i    
    helper = SwaptionHelper(
        Period(t, Years), 
        Period(tenor, Years), 
        QuoteHandle(SimpleQuote(vol[i])), 
        USDLibor(Period(3, Months), curveHandle), 
        Period(1, Years), 
        dayCounter, 
        dayCounter, 
        curveHandle)    
    calibrator.AddCalibrationHelper(helper)

# create model and pricing engine, calibrate model and print calibrated parameters
print('case 1 : calibrate all involved parameters (HW1F : reversion, sigma)')
model = HullWhite(curveHandle)
engine = JamshidianSwaptionEngine(model)
fixedParameters = []
calibrator.Calibrate(model, engine, curveHandle, fixedParameters)
print('calibrated reversion: ' + str(round(model.params()[0], 5)))
print('calibrated sigma: ' + str(round(model.params()[1], 5)))
print()

print('case 2 : calibrate sigma and fix reversion to 0.05')
model = HullWhite(curveHandle, 0.05, 0.0001)
engine = JamshidianSwaptionEngine(model)
fixedParameters = [True, False]
calibrator.Calibrate(model, engine, curveHandle, fixedParameters)
print('fixed reversion: ' + str(round(model.params()[0], 5)))
print('calibrated sigma: ' + str(round(model.params()[1], 5)))
print()

print('case 3 : calibrate reversion and fix sigma to 0.01')
model = HullWhite(curveHandle, 0.05, 0.01)
engine = JamshidianSwaptionEngine(model)
fixedParameters = [False, True]
calibrator.Calibrate(model, engine, curveHandle, fixedParameters)
print('calibrated reversion: ' + str(round(model.params()[0], 5)))
print('fixed sigma: ' + str(round(model.params()[1], 5)))
print()

QuantLib : Hull-White one-factor model calibration

Sometimes during the last year I published one post on simulating Hull-White interest rate paths using Quantlib. My conclusion was, that with all the tools provided by this wonderful library, this task should be (relatively) easy thing to do. However, as we know the devil is always in the details - in this case, in the actual process parameters (mean reversion, sigma). Before processing any simulation, we have to get those parameters from somewhere. In this post, we use Quantlib tools for calibrating our model parameters to swaption prices existing in the market.

There are so many variations of this model out there (depending on time-dependencies of different parameters), but the following is the model we are going to use in this example.

Parameters alpha and sigma are constants and theta is time-dependent variable (usually calibrated to current yield curve). Swaption surface is presented in the picture below. Within the table below, time to swaption maturity has been defined in vertical axis, while tenor for underlying swap contract has been defined in horizontal axis. Co-terminal swaptions (used in calibration process) have been specifically marked with yellow colour.

Calibration results

Test cases for three different calibration schemes are included in the example program. More specifically, we :

1. calibrate the both parameters
2. calibrate sigma parameter and freeze reversion parameter to 0.05
3. calibrate reversion parameter and freeze sigma parameter to 0.01

With the given data, we get the following results for these calibration schemes.

The program

As preparatory task, builder class for constructing yield curve should be implemented into a new project from here. After this task, the following header file (ModelCalibrator.h) and tester file (Tester.cpp) should be added to this new project.

First, piecewise yield curve (USD swap curve) and swaption volatilities (co-terminal swaptions) are created by using two free functions in tester. All the data has been hard-coded inside these functions. Needless to say, in any production-level program, this data feed should come from somewhere else. After required market data and ModelCalibrator object has been created, calibration helpers for diagonal swaptions are going to be created and added to ModelCalibrator. Finally, test cases for different calibration schemes are processed.

// ModelCalibrator.h
#pragma once
#include <ql\quantlib.hpp>
#include <algorithm>
//
namespace MJModelCalibratorNamespace
{
 using namespace QuantLib;
 //
 template <typename MODEL, typename OPTIMIZER = LevenbergMarquardt>
 class ModelCalibrator
 {
 public:
  // default implementation (given values) for end criteria class
  ModelCalibrator(const EndCriteria& endCriteria = EndCriteria(1000, 500, 0.0000001, 0.0000001, 0.0000001))
   : endCriteria(endCriteria)
  { }
  void AddCalibrationHelper(boost::shared_ptr<CalibrationHelper>& helper)
  {
   // add any type of calibration helper
   helpers.push_back(helper);
  }
  void Calibrate(boost::shared_ptr<MODEL>& model,
   const boost::shared_ptr<PricingEngine>& pricingEngine,
   const Handle<YieldTermStructure>& curve,
   const std::vector<bool> fixedParameters = std::vector<bool>())
  {
   // assign pricing engine to all calibration helpers
   std::for_each(helpers.begin(), helpers.end(),
    [&pricingEngine](boost::shared_ptr<CalibrationHelper>& helper)
   { helper->setPricingEngine(pricingEngine); });
   //
   // create optimization model for calibrating requested model
   OPTIMIZER solver;
   //
   if (fixedParameters.empty())
   {
    // calibrate all involved parameters
    model->calibrate(helpers, solver, this->endCriteria);
   }
   else
   {
    // calibrate all involved non-fixed parameters
    // hard-coded : vector for weights and constraint type
    model->calibrate(helpers, solver, this->endCriteria, NoConstraint(), std::vector<Real>(), fixedParameters);
   }
  }
 private:
  EndCriteria endCriteria;
  std::vector<boost::shared_ptr<CalibrationHelper>> helpers;
 };
}
//
//
//
//
// Tester.cpp
#include "PiecewiseCurveBuilder.cpp"
#include "ModelCalibrator.h"
#include <iostream>
//
using namespace MJPiecewiseCurveBuilderNamespace;
namespace MJ_Calibrator = MJModelCalibratorNamespace;
//
// declarations for two free functions used to construct required market data
void CreateCoTerminalSwaptions(std::vector<Volatility>& diagonal);
void CreateYieldCurve(RelinkableHandle<YieldTermStructure>& curve, 
 const Date& settlementDate, const Calendar& calendar);
//
//
//
int main()
{
 try
 {
  // dates
  Date tradeDate(4, September, 2017);
  Settings::instance().evaluationDate() = tradeDate;
  Calendar calendar = TARGET();
  Date settlementDate = calendar.advance(tradeDate, Period(2, Days), ModifiedFollowing);
  //
  // market data : create piecewise yield curve
  RelinkableHandle<YieldTermStructure> curve;
  CreateYieldCurve(curve, settlementDate, calendar);
  //
  // market data : create co-terminal swaption volatilities
  std::vector<Volatility> diagonal;
  CreateCoTerminalSwaptions(diagonal);
  //
  // create model calibrator
  MJ_Calibrator::ModelCalibrator<HullWhite> modelCalibrator;
  //
  // create and add calibration helpers to model calibrator
  boost::shared_ptr<IborIndex> floatingIndex(new USDLibor(Period(3, Months), curve));
  for (unsigned int i = 0; i != diagonal.size(); ++i)
  {
   int timeToMaturity = i + 1;
   int underlyingTenor = diagonal.size() - i;
   //
   // using 1st constructor for swaption helper class
   modelCalibrator.AddCalibrationHelper(boost::shared_ptr<CalibrationHelper>(new SwaptionHelper(
     Period(timeToMaturity, Years), // time to swaption maturity
     Period(underlyingTenor, Years), // tenor of the underlying swap
     Handle<Quote>(boost::shared_ptr<Quote>(new SimpleQuote(diagonal[i]))), // swaption volatility
     floatingIndex, // underlying floating index
     Period(1, Years), // tenor for underlying fixed leg
     Actual360(), // day counter for underlying fixed leg
     floatingIndex->dayCounter(), // day counter for underlying floating leg
     curve))); // term structure
  }
  //
  // create model and pricing engine, calibrate model and print calibrated parameters
  // case 1 : calibrate all involved parameters (HW1F : reversion, sigma)
  boost::shared_ptr<HullWhite> model(new HullWhite(curve));
  boost::shared_ptr<PricingEngine> jamshidian(new JamshidianSwaptionEngine(model));
  modelCalibrator.Calibrate(model, jamshidian, curve);
  std::cout << "calibrated reversion = " << model->params()[0] << std::endl;
  std::cout << "calibrated sigma = " << model->params()[1] << std::endl;
  std::cout << std::endl;
  //
  // case 2 : calibrate sigma and fix reversion to famous 0.05
  model = boost::shared_ptr<HullWhite>(new HullWhite(curve, 0.05, 0.0001));
  jamshidian = boost::shared_ptr<PricingEngine>(new JamshidianSwaptionEngine(model));
  std::vector<bool> fixedReversion = { true, false };
  modelCalibrator.Calibrate(model, jamshidian, curve, fixedReversion);
  std::cout << "fixed reversion = " << model->params()[0] << std::endl;
  std::cout << "calibrated sigma = " << model->params()[1] << std::endl;
  std::cout << std::endl;
  //
  // case 3 : calibrate reversion and fix sigma to 0.01
  model = boost::shared_ptr<HullWhite>(new HullWhite(curve, 0.05, 0.01));
  jamshidian = boost::shared_ptr<PricingEngine>(new JamshidianSwaptionEngine(model));
  std::vector<bool> fixedSigma = { false, true };
  modelCalibrator.Calibrate(model, jamshidian, curve, fixedSigma);
  std::cout << "calibrated reversion = " << model->params()[0] << std::endl;
  std::cout << "fixed sigma = " << model->params()[1] << std::endl;
 }
 catch (std::exception& e)
 {
  std::cout << e.what() << std::endl;
 }
 return 0;
}
//
void CreateCoTerminalSwaptions(std::vector<Volatility>& diagonal)
{
 // hard-coded data
 // create co-terminal swaptions 
 diagonal.push_back(0.3133); // 1x10
 diagonal.push_back(0.3209); // 2x9
 diagonal.push_back(0.3326); // 3x8
 diagonal.push_back(0.331); // 4x7
 diagonal.push_back(0.3281); // 5x6
 diagonal.push_back(0.318); // 6x5
 diagonal.push_back(0.3168); // 7x4
 diagonal.push_back(0.3053); // 8x3
 diagonal.push_back(0.2992); // 9x2
 diagonal.push_back(0.3073); // 10x1
}
//
void CreateYieldCurve(RelinkableHandle<YieldTermStructure>& curve,
 const Date& settlementDate, const Calendar& calendar) 
{
 // hard-coded data
 // create piecewise yield curve by using builder class
 DayCounter curveDaycounter = Actual360();
 PiecewiseCurveBuilder<ZeroYield, Linear> builder;
 //
 // cash rates
 pQuote q_1W(new SimpleQuote(0.012832));
 pIndex i_1W(new USDLibor(Period(1, Weeks)));
 builder.AddDeposit(q_1W, i_1W);
 //
 pQuote q_1M(new SimpleQuote(0.012907));
 pIndex i_1M(new USDLibor(Period(1, Months)));
 builder.AddDeposit(q_1M, i_1M);
 //
 pQuote q_3M(new SimpleQuote(0.0131611));
 pIndex i_3M(new USDLibor(Period(3, Months))); 
 builder.AddDeposit(q_3M, i_3M);
 //
 // futures
 Date IMMDate;
 pQuote q_DEC17(new SimpleQuote(98.5825));
 IMMDate = IMM::nextDate(settlementDate + Period(3, Months));
 builder.AddFuture(q_DEC17, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_MAR18(new SimpleQuote(98.5425));
 IMMDate = IMM::nextDate(settlementDate + Period(6, Months));
 builder.AddFuture(q_MAR18, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_JUN18(new SimpleQuote(98.4975));
 IMMDate = IMM::nextDate(settlementDate + Period(9, Months));
 builder.AddFuture(q_JUN18, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_SEP18(new SimpleQuote(98.4475));
 IMMDate = IMM::nextDate(settlementDate + Period(12, Months));
 builder.AddFuture(q_SEP18, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_DEC18(new SimpleQuote(98.375));
 IMMDate = IMM::nextDate(settlementDate + Period(15, Months));
 builder.AddFuture(q_DEC18, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_MAR19(new SimpleQuote(98.3425));
 IMMDate = IMM::nextDate(settlementDate + Period(18, Months));
 builder.AddFuture(q_MAR19, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_JUN19(new SimpleQuote(98.3025));
 IMMDate = IMM::nextDate(settlementDate + Period(21, Months));
 builder.AddFuture(q_JUN19, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_SEP19(new SimpleQuote(98.2675));
 IMMDate = IMM::nextDate(settlementDate + Period(24, Months));
 builder.AddFuture(q_SEP19, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_DEC19(new SimpleQuote(98.2125));
 IMMDate = IMM::nextDate(settlementDate + Period(27, Months));
 builder.AddFuture(q_DEC19, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_MAR20(new SimpleQuote(98.1775));
 IMMDate = IMM::nextDate(settlementDate + Period(30, Months));
 builder.AddFuture(q_MAR20, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_JUN20(new SimpleQuote(98.1425));
 IMMDate = IMM::nextDate(settlementDate + Period(33, Months));
 builder.AddFuture(q_JUN20, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 // swaps
 pIndex swapFloatIndex(new USDLibor(Period(3, Months))); 
 pQuote q_4Y(new SimpleQuote(0.01706));
 builder.AddSwap(q_4Y, Period(4, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_5Y(new SimpleQuote(0.0176325));
 builder.AddSwap(q_5Y, Period(5, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_6Y(new SimpleQuote(0.01874));
 builder.AddSwap(q_6Y, Period(6, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_7Y(new SimpleQuote(0.0190935));
 builder.AddSwap(q_7Y, Period(7, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_8Y(new SimpleQuote(0.02011));
 builder.AddSwap(q_8Y, Period(8, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_9Y(new SimpleQuote(0.02066));
 builder.AddSwap(q_9Y, Period(9, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_10Y(new SimpleQuote(0.020831));
 builder.AddSwap(q_10Y, Period(10, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_11Y(new SimpleQuote(0.02162));
 builder.AddSwap(q_11Y, Period(11, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_12Y(new SimpleQuote(0.0217435));
 builder.AddSwap(q_12Y, Period(12, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_15Y(new SimpleQuote(0.022659));
 builder.AddSwap(q_15Y, Period(15, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_20Y(new SimpleQuote(0.0238125));
 builder.AddSwap(q_20Y, Period(20, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_25Y(new SimpleQuote(0.0239385));
 builder.AddSwap(q_25Y, Period(25, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_30Y(new SimpleQuote(0.02435));
 builder.AddSwap(q_30Y, Period(30, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 curve = builder.GetCurveHandle(settlementDate, curveDaycounter);
}

I have noticed, that there are some rule-of-thumbs, but the actual calibration for any of these models is not so straightforward as one may think. There are some subtle issues on market data quality and products under pricing, which are leaving us with relatively high degree of freedom. Further step towards the calibration red pill can be taken by checking out this excellent research paper, written by the guys from Mizuho Securities.

Finally, as always, thanks a lot again for reading this blog.
-Mike

AlgLib : Ho-Lee Calibration Using Levenberg-Marquardt algorithm in VBA

Some time ago, I published one possible C# implementation for Ho-Lee one-factor model calibration scheme using AlgLib numerical libraries. This time I will present an implementation for the same scheme in VBA. General information concerning Levenberg-Marquardt algorithm implementation in AlgLib has been presented in here. Libraries for VBA (collection of BAS module files, which are going to be included into VBA project) can be downloaded from here.

A few words about this implementation. AlgLibLMASolver class uses AlgLib library functions (functions from 21 different modules) for processing (creating optimization model, setting conditions, processing iterations). One data member within this class is having a type of IModel. This data member is actually a reference to an interface, which provides a set of functions for all required calculations (objective function value, values for function terms, partial derivative values for function terms). Since all possible implementations for any interface method must honor signatures exactly, there is a problem with VBA since it does not have a real constructor mechanism. I have chewed this issue in here. It might help to explain the reason, why I have been distributing input parameters for interface implementation by using Dictionary object. Finally, HoLeeZeroCouponCalibration class is implementing IModel interface (a set of functions for all required calculations). In essence, algorithms (AlgLib-related processing) and data (values calculated specifically by using Ho-Lee model) have been completely separated. Needless to say, this type of scheme is flexible for new implementations.

Create a new VBA project and copy-paste the following classes and modules into this project. Also, import all required 21 AlgLib BAS files into this project.

' CLASS : AlgLibLMASolver
Option Explicit
'
' The following 21 AlgLib modules are required for succesfull compilation of this project :
' ablas, ablasf, ap, bdsvd, blas, creflections, densesolver, hblas, linmin, matinv,
' minlbfgs, minlm, ortfac, rcond, reflections, rotations, safesolve, sblas, svd, trfac, xblas
'
Private state As MinLMState
Private report As MinLMReport
Private n As Long
Private m As Long
Private x() As Double
Private model As IModel
Private epsF As Double
Private epsG As Double
Private epsX As Double
Private iterations As Long
'
Public Function initialize( _
    ByVal numberOfVariables As Long, _
    ByVal numberOfEquations As Long, _
    ByRef changingVariables() As Double, _
    ByRef callbackModel As IModel, _
    ByVal epsilonF As Double, _
    ByVal epsilonG As Double, _
    ByVal epsilonX As Double, _
    ByVal maximumIterations As Long)
    '
    n = numberOfVariables
    m = numberOfEquations
    x = changingVariables
    Set model = callbackModel
    epsF = epsilonF
    epsG = epsilonG
    epsX = epsilonX
    iterations = maximumIterations
End Function
'
Public Sub Solve()
    '
    ' create solver scheme using functions and analytical partial derivatives
    Call MinLMCreateFJ(n, m, x, state)
    ' set stopping conditions
    Call MinLMSetCond(state, epsG, epsF, epsX, iterations)
    '
    ' process iterations
    Do While MinLMIteration(state)
        '
        ' calculate value for objective function
        If (state.NeedF) Then
            '
            model.callBackObjectiveFunction state
        End If
        '
        ' calculate values for functions and partial derivatives
        If (state.NeedFiJ) Then
            '
            model.callBackFunction state
            model.callBackJacobian state
        End If
    Loop
    '
    ' process results
    Call MinLMResults(state, x, report)
End Sub
'
' public accessor to (MinLMState) state
Public Property Get GetState() As MinLMState
    GetState = state
End Property
'
' public accessor to (MinLMReport) report
Public Property Get GetReport() As MinLMReport
    GetReport = report
End Property
'
' public accessor to hard-coded report
Public Property Get GetPrettyPrintReport() As String
    '
    Dim message As String
    message = "*** AlgLibLMASolver execution report " + VBA.CStr(VBA.Now) + " ***" + VBA.vbNewLine
    message = message + "TerminationType : " + VBA.CStr(report.TerminationType) + VBA.vbNewLine
    message = message + "Iterations : " + VBA.CStr(report.IterationsCount) + VBA.vbNewLine
    message = message + "Objective function : " + VBA.CStr(state.f) + VBA.vbNewLine
    message = message + VBA.vbNewLine
    '
    Dim i As Integer
    For i = 0 To (state.n - 1)
        message = message + "x(" + VBA.CStr(i) + ") = " + VBA.CStr(state.x(i)) + VBA.vbNewLine
    Next i
    '
    GetPrettyPrintReport = message
End Property
'
'
'
'
'
' CLASS : IModel
Option Explicit
' set of functions for IModel interface
Public Function initialize(ByRef parameters As Scripting.Dictionary)
    ' assign required member data wrapped into dictionary
End Function
'
Public Function callBackObjectiveFunction(ByRef state As MinLMState)
    ' calculate objective function value
End Function
'
Public Function callBackFunction(ByRef state As MinLMState)
    ' calculate values for (non-squared) function terms
End Function
'
Public Function callBackJacobian(ByRef state As MinLMState)
    ' calculate partial derivative values for (non-squared) function terms
End Function
'
'
'
'
'
' CLASS : HoLeeZeroCouponCalibration
Option Explicit
'
Implements IModel
'
Private s As Double
Private r As Double
Private t() As Double
Private z() As Double
'
Private Function IModel_initialize(ByRef parameters As Scripting.IDictionary)
    '
    s = parameters(HOLEE_PARAMETERS.sigma)
    r = parameters(HOLEE_PARAMETERS.shortRate)
    t = parameters(HOLEE_PARAMETERS.maturity)
    z = parameters(HOLEE_PARAMETERS.zeroCouponBond)
End Function
'
Private Function IModel_callBackObjectiveFunction(ByRef state As MinLMState)
    '
    ' calculate value for aggregate objective function
    Dim i As Integer
    Dim hoLeeZero As Double
    Dim f As Double: f = 0
    '
    ' loop through number of equations
    For i = 0 To (state.m - 1)
        '
        hoLeeZero = VBA.Exp(-(1 / 2) * state.x(i) * (t(i) ^ 2) + (1 / 6) * (s ^ 2) * (t(i) ^ 3) - r * t(i))
        f = f + (z(i) - hoLeeZero) ^ 2
    Next i
    state.f = f
End Function
'
Private Function IModel_callBackFunction(ByRef state As MinLMState)
    '
    ' calculate values for (non-squared) function terms
    Dim i As Integer
    Dim hoLeeZero As Double
    '
    ' loop through number of equations
    For i = 0 To (state.m - 1)
        '
        hoLeeZero = VBA.Exp(-(1 / 2) * state.x(i) * (t(i) ^ 2) + (1 / 6) * (s ^ 2) * (t(i) ^ 3) - r * t(i))
        state.FI(i) = (z(i) - hoLeeZero)
    Next i
End Function
'
Private Function IModel_callBackJacobian(ByRef state As MinLMState)
    '
    ' calculate partial derivative values for (non-squared) function terms
    Dim i As Integer, J As Integer
    Dim hoLeeZero As Double
    '
    ' 1. individual (non-squared) function terms
    ' loop through number of equations
    For i = 0 To (state.m - 1)
        '
        hoLeeZero = VBA.Exp(-(1 / 2) * state.x(i) * (t(i) ^ 2) + (1 / 6) * (s ^ 2) * (t(i) ^ 3) - r * t(i))
        state.FI(i) = (z(i) - hoLeeZero)
    Next i
    '
    ' 2. partial derivatives for all (non-squared) function terms
    ' loop through number of equations
    For i = 0 To (state.m - 1)
        '
    ' loop through number of variables
        For J = 0 To (state.n - 1)
            '
            Dim derivative As Double: derivative = 0
            ' partial derivative is non-zero only for diagonal cases
            If (i = J) Then
                derivative = (1 / 2) * VBA.Exp(1) * t(J) ^ 2
                state.J(i, J) = derivative
            End If
        Next J
    Next i
End Function
'
'
'
'
'
' MODULE : DataStructures
Option Explicit
'
Public Enum HOLEE_PARAMETERS
    sigma = 1
    shortRate = 2
    maturity = 3
    zeroCouponBond = 4
End Enum
'
'
'
'
'
' TESTER MODULE
Option Explicit
'
' Ho-Lee model calibration example
Public Sub AlglibTester()
    '
    ' MODEL part
    ' construct all required inputs and model to be calibrated
    Dim sigma As Double: sigma = 0.00039
    Dim shortRate As Double: shortRate = 0.00154
    '
    Dim maturity(0 To 9) As Double
    maturity(0) = 1: maturity(1) = 2: maturity(2) = 3: maturity(3) = 4: maturity(4) = 5:
    maturity(5) = 6: maturity(6) = 7: maturity(7) = 8: maturity(8) = 9: maturity(9) = 10
    '
    Dim zero(0 To 9) As Double
    zero(0) = 0.9964: zero(1) = 0.9838: zero(2) = 0.9611: zero(3) = 0.9344: zero(4) = 0.9059:
    zero(5) = 0.8769: zero(6) = 0.8478: zero(7) = 0.8189: zero(8) = 0.7905: zero(9) = 0.7626
    '
    ' assign parameters into dictionary wrapper
    Dim parameters As New Scripting.Dictionary
    parameters.Add HOLEE_PARAMETERS.sigma, sigma
    parameters.Add HOLEE_PARAMETERS.shortRate, shortRate
    parameters.Add HOLEE_PARAMETERS.maturity, maturity
    parameters.Add HOLEE_PARAMETERS.zeroCouponBond, zero
    '
    ' create and initialize calibration model
    Dim model As IModel: Set model = New HoLeeZeroCouponCalibration
    model.initialize parameters
    '
    ' SOLVER part
    Dim Theta(0 To 9) As Double ' assign initial guesses
    Theta(0) = 0.001: Theta(1) = 0.001: Theta(2) = 0.001: Theta(3) = 0.001: Theta(4) = 0.001:
    Theta(5) = 0.001: Theta(6) = 0.001: Theta(7) = 0.001: Theta(8) = 0.001: Theta(9) = 0.001
    '
    Dim numberOfVariables As Long: numberOfVariables = 10
    Dim numberOfEquations As Long: numberOfEquations = 10
    Dim epsilonF As Double: epsilonF = 0.000000000001
    Dim epsilonG As Double: epsilonG = 0.000000000001
    Dim epsilonX As Double: epsilonX = 0.000000000001
    Dim maximumIterations As Long: maximumIterations = 25000
    '
    ' create and initialize solver model
    Dim solver As New AlgLibLMASolver
    solver.initialize _
        numberOfVariables, _
        numberOfEquations, _
        Theta, _
        model, _
        epsilonF, _
        epsilonG, _
        epsilonX, _
        maximumIterations
    '
    ' solve calibration model
    solver.Solve
    '
    ' print hard-coded report containing values for
    ' objective function, variables and other information
    Debug.Print solver.GetPrettyPrintReport
End Sub
'

The results from this calibration model have been verified against the previous results.

Importing several files into project may involve considerable amount of cruel and unusual repetitive labour. For this specific reason, I have also been woodshedding a separate module (employing VBIDE object), which might give some relief when babysitting those AlgLib modules.

' The following dll libraries need to be referenced :
' Microsoft Visual Basic for Applications Extensibility 5.X, Microsoft Scripting Runtime
Option Explicit
Option Base 0
'
' address to a list, which contains all BAS files which will be included into project
Const listFolderPathName As String = "C:\AlgLib\vba\AlgLibLMAModules.txt"
' address to a folder, which contains all AlgLib BAS files
Const moduleFolderPathName  As String = "C:\AlgLib\vba\alglib-2.6.0.vb6\vb6\src\"
' select TRUE, if Require Variable Declaration in editor is tagged
Const removeOptionExplicitDirective As Boolean = True
'
Public Sub ImportModules()
    '
    ' create a list of modules to be imported into this project
    Dim list() As String: list = createProjectModuleList
    ' import modules into active project
    import list
End Sub
'
Public Sub ExportModules()
    '
    ' create a list of modules to be exported from this project
    Dim list() As String: list = createProjectModuleList
    ' export modules from active project into a defined folder
    export list
End Sub
'
Public Sub RemoveModules()
    '
    ' create a list of modules to be removed from this project
    Dim list() As String: list = createProjectModuleList
    ' delete modules from active project
    remove list
End Sub
'
Private Function import(ByRef list() As String)
    '
    Dim editor As VBIDE.VBProject
    Set editor = ActiveWorkbook.VBProject
    Dim fileSystem As Scripting.FileSystemObject: Set fileSystem = New Scripting.FileSystemObject
    '
    ' loop through all files in a specific source folder for modules
    Dim filesInGivenList As Integer: filesInGivenList = UBound(list) + 1
    If (filesInGivenList = 0) Then Exit Function
    '
    Dim module As VBIDE.VBComponent
    Dim file As Scripting.file
    For Each file In fileSystem.GetFolder(moduleFolderPathName).Files
        '
        ' if there is a given list of specific files to be included
        ' select only the files in that list to be imported into project
        If Not (moduleIsIncluded(file.Name, list)) Then GoTo skipPoint
        '
        Set module = editor.VBComponents.Add(vbext_ct_StdModule)
        If (removeOptionExplicitDirective) Then module.CodeModule.DeleteLines 1
        module.Name = VBA.Split(file.Name, ".")(0)
        module.CodeModule.AddFromFile file.Path
skipPoint:
    Next
End Function
'
Private Function export(ByRef list() As String)
    '
    Dim filesInGivenList As Integer: filesInGivenList = UBound(list) + 1
    If (filesInGivenList = 0) Then Exit Function
    '
    Dim editor As VBIDE.VBProject
    Set editor = ActiveWorkbook.VBProject
    Dim module As VBIDE.VBComponent
    '
    ' loop through all modules
    For Each module In editor.VBComponents
        '
        ' export module only if it is included in the list
        If (moduleIsIncluded(module.Name + ".bas", list)) Then
            module.export moduleFolderPathName + module.Name + ".bas"
        End If
    Next
End Function
'
Private Function remove(ByRef list() As String)
    '
    Dim filesInGivenList As Integer: filesInGivenList = UBound(list) + 1
    If (filesInGivenList = 0) Then Exit Function
    '
    Dim editor As VBIDE.VBProject
    Set editor = ActiveWorkbook.VBProject
    Dim module As VBIDE.VBComponent
    '
    ' loop through all modules
    For Each module In editor.VBComponents
        '
        ' remove module only if it is included in the list
        If (moduleIsIncluded(module.Name + ".bas", list)) Then
            module.Collection.remove module
        End If
    Next
End Function
'
Private Function moduleIsIncluded(ByVal FileName As String, ByRef list() As String) As Boolean
    '
    ' check if a given file name is in the list
    Dim isIncluded As Boolean: isIncluded = False
    Dim i As Integer
    For i = 0 To UBound(list)
        If (FileName = list(i)) Then
            isIncluded = True
            Exit For
        End If
    Next i
    moduleIsIncluded = isIncluded
End Function
'
Private Function createProjectModuleList() As String()
    '
    ' create a list of file names from text file
    Dim fileSystem As Scripting.FileSystemObject: Set fileSystem = New Scripting.FileSystemObject
    Dim fileReader As Scripting.TextStream: Set fileReader = fileSystem.OpenTextFile(listFolderPathName, ForReading)
    Dim fileStreams As String: fileStreams = fileReader.ReadAll
    Dim streams As Variant: streams = VBA.Split(fileStreams, VBA.vbNewLine)
    Dim list() As String: ReDim list(0 To UBound(streams))
    Dim i As Integer
    For i = 0 To UBound(streams)
        list(i) = VBA.Trim(streams(i))
    Next i
    createProjectModuleList = list
End Function
'

Finally, thanks for reading this blog.
-Mike

Alglib : SABR calibration in C#

Discovering new tools for implementing numerical stuff is refreshing. This time, I wanted to present my calibration implementation for SABR model in C# by using Alglib numerical library. The thing I really like in Alglib, is its easiness of use : just find Alglib download page, download proper C# version, create a new C# project and add reference to DLL file included in download folder - that's all. Like in my previous post, optimization routine for this task has also been implemented by using Levenberg-Marquardt algorithm.

FLYING CARPETS

Let us see the end product first. SABR coefficients have been calibrated for 16 different volatility skews. By using four calibrated parameters within SABR model, we should get extremely close replica of original volatility skew. Original forward volatility surface and its SABR-calibrated replica is presented in the screenshot below. Source data has been taken from this great book (Excel files included) written by Richard Flavell. The essence of SABR model calibration has been chewed in the chapter 10.

SKEW DISSECTION

In order to validate this C# implementation, I have been replicating example (Non-shifted SABR calibration results) found in MathWorks page on SABR model calibration. Original example data and calibration results of this C# implementation are presented in the screenshot below.

Additional results have been presented for shifted SABR model calibrations. On the current low rate environment (in some currencies) we may enter into negative strike/forward rate territory. As this happens, log-normal model will crash. For such cases, using shifted log-normal SABR model is deceptively simple : constant shift will be added into strike prices and forward rates, in order to be able to calibrate the model.

It should be noted, that (in this implementation) beta coefficient is assumed to be a given estimated constant as this is usually a market practice. Moreover, the program is able to calibrate the model by solving three remaining coefficients (alpha, rho, nu) directly or then solving alpha implicitly from two remaining coefficients (rho, nu). These two calibration schemes have been presented in MathWorks example (Method 1 and Method 2). Finally, calibration speed (this should not be a problem anyway) and especially accuracy can be greatly improved by using Vector-Jacobian scheme (done in my previous post) and explicitly defined analytical partial derivatives instead of their numerical equivalents (finite difference).

THE PROGRAM

In order to test this calibration program, create a new C# project (SABRCalibrationTester), copy-paste the following code into a new CS file and add reference to Alglib DLL file.

using System;
using System.Collections.Generic;
using System.Linq;

namespace SABRCalibrationTester
{
    class Program
    {
        static void Main(string[] args)
        {
            double beta = 0.5;
            double expiry = 3.0027397260274;
            double forward = 0.035;
            double atmVolatility = 0.366;
            double[] strike = new double[] { 0.02, 0.025, 0.03, 0.035, 0.04, 0.045, 0.05 };
            double[] volatility = new double[] { 0.456, 0.416, 0.379, 0.366, 0.378, 0.392, 0.4 };
            double[] coefficients = null;
            //
            // case 1 : use explicit alpha
            // order of initial coefficients : alpha, rho, nu
            SABR skew_3Y = new SABR(beta, expiry, forward, atmVolatility, strike, volatility);
            coefficients = new double[] { 0.25, 0.25, 0.25 };
            LevenbergMarquardtSolver.Solve(ref coefficients, 7, (IFunctionVector)skew_3Y);
            Print(skew_3Y, "explicit alpha, no shift");
            //
            // case 2 : use implicit alpha
            // order of initial coefficients : rho, nu
            coefficients = new double[] { 0.25, 0.25 };
            LevenbergMarquardtSolver.Solve(ref coefficients, 7, (IFunctionVector)skew_3Y);
            Print(skew_3Y, "implicit alpha, no shift");
            //
            // case 3 : use explicit alpha and 4% shift factor
            // order of initial coefficients : alpha, rho, nu
            // create negative forward/strike scenario by shifting forward rate and strike prices down by 4%
            forward = -0.005;
            strike = new double[] { -0.02, -0.015, -0.01, -0.005, 0.0, 0.005, 0.01 };
            double shift = 0.04;
            SABR shiftedSkew_3Y = new SABR(beta, expiry, forward, atmVolatility, strike, volatility, shift);
            coefficients = new double[] { 0.25, 0.25, 0.25 };
            LevenbergMarquardtSolver.Solve(ref coefficients, 7, (IFunctionVector)shiftedSkew_3Y);
            Print(shiftedSkew_3Y, "explicit alpha, 4% shift");
            //
            // case 4 : use implicit alpha and 4% shift factor
            // order of initial coefficients : rho, nu
            coefficients = new double[] { 0.25, 0.25 };
            LevenbergMarquardtSolver.Solve(ref coefficients, 7, (IFunctionVector)shiftedSkew_3Y);
            Print(shiftedSkew_3Y, "implicit alpha, 4% shift");
        }
        public static void Print(SABR sabr, string message)
        {
            // print message, SABR-calibrated coefficients and volatility skew
            Console.WriteLine(message);
            Console.WriteLine(String.Format("Alpha = {0:0.000%}", sabr.Alpha));
            Console.WriteLine(String.Format("Beta = {0:0.000%}", sabr.Beta));
            Console.WriteLine(String.Format("Rho = {0:0.000%}", sabr.Rho));
            Console.WriteLine(String.Format("Nu = {0:0.000%}", sabr.Nu));
            List<double> skew = sabr.SABRSkew().ToList<double>();
            skew.ForEach(it => Console.WriteLine(String.Format("{0:0.000%}", it)));
            Console.WriteLine();
        }
    }
    //
    //
    //
    //
    /// <summary>
    /// interface for function vector callback method, required by LevenbergMarquardtSolver class
    /// </summary>
    public interface IFunctionVector
    {
        void callback(double[] x, double[] fi, object obj);
    }
    //
    /// <summary>
    /// interface : definition for jacobian matrix callback method, optionally required by LevenbergMarquardtSolver class
    /// </summary>
    public interface IJacobianMatrix
    {
        void callback(double[] x, double[] fi, double[,] jac, object obj);
    }
    //
    /// <summary>
    /// static wrapper class for Alglib Levenberg-Marquardt solver implementation
    /// </summary>
    public static class LevenbergMarquardtSolver
    {
        // alglib report is exposed to public
        public static alglib.minlmreport report;
        private static alglib.minlmstate state;
        //
        /// <summary>
        /// create LMA minimization solver using only function vector ('V' vector method)
        /// </summary>
        /// <param name="x">array of initial coefficients</param>
        /// <param name="m">number of functions</param>
        /// <param name="functionVector">IFunctionVector interface implementation</param>
        public static void Solve(ref double[] x, int m, IFunctionVector functionVector,
            // optional parameters which are set to default values
            double diffstep = 1E-08, double epsg = 1E-15, double epsf = 0.0, double epsx = 0.0, int maxits = 0)
        {
            // create LM model, set termination conditions, perform optimization and get results
            alglib.minlmcreatev(m, x, diffstep, out state);
            alglib.minlmsetcond(state, epsg, epsf, epsx, maxits);
            alglib.minlmoptimize(state, functionVector.callback, null, null);
            alglib.minlmresults(state, out x, out report);
        }
        //
        /// <summary>
        /// create LMA minimization solver using function vector and Jacobian matrix ('VJ' vector-jacobian method)
        /// </summary>
        /// <param name="x">array of initial coefficients</param>
        /// <param name="m">number of functions</param>
        /// <param name="functionVector">IFunctionVector interface implementation</param>
        /// <param name="jacobianMatrix">IJacobianMatrix interface implementation</param>
        public static void Solve(ref double[] x, int m, IFunctionVector functionVector, IJacobianMatrix jacobianMatrix,
            // optional parameters which are set to default values
            double epsg = 1E-15, double epsf = 0.0, double epsx = 0.0, int maxits = 0)
        {
            // create LM model, set termination conditions, perform optimization and get results
            alglib.minlmcreatevj(m, x, out state);
            alglib.minlmsetcond(state, epsg, epsf, epsx, maxits);
            alglib.minlmoptimize(state, functionVector.callback, jacobianMatrix.callback, null, null);
            alglib.minlmresults(state, out x, out report);
        }
    }
    //
    //
    //
    //
    /// <summary>
    /// class calibrating SABR volatility skew for a given set of parameters
    /// class implements IFunctionVector interface, uses static MathTools class
    /// </summary>
    public class SABR : IFunctionVector
    {
        // private member data
        private double alpha;
        private double beta;
        private double rho;
        private double nu;
        private double expiry;
        private double forward;
        private double atmVolatility;
        private double shift;
        private double[] strike;
        private double[] marketVolatility;
        //
        // public read-only access to SABR coefficients
        public double Alpha { get { return this.alpha; } }
        public double Beta { get { return this.beta; } }
        public double Rho { get { return this.rho; } }
        public double Nu { get { return this.nu; } }
        public double Expiry { get { return this.expiry; } }
        public double Forward { get { return this.forward; } }
        public double ATMVolatility { get { return this.atmVolatility; } }
        public double Shift { get { return this.shift; } }
        public double[] Strike { get { return this.strike; } }
        public double[] MarketVolatility { get { return this.marketVolatility; } }
        //
        /// <summary>
        /// constructor
        /// </summary>
        /// <param name="beta">assumed to be a given constant within interval [0, 1]</param>
        /// <param name="expiry">expiration time for forward</param>
        /// <param name="forward">forward market quote</param>
        /// <param name="atmVolatility">at-the-money market volatility quote</param>
        /// <param name="strike">array of strike prices for corresponding market volatility quotes</param>
        /// <param name="volatility">array of market volatility quotes</param>
        public SABR(double beta, double expiry, double forward, double atmVolatility, double[] strike, double[] volatility, double shift = 0.0)
        {
            this.beta = beta;
            this.expiry = expiry;
            this.forward = forward;
            this.atmVolatility = atmVolatility;
            this.strike = strike;
            this.marketVolatility = volatility;
            this.shift = shift;
        }
        /// <summary>
        /// retrieve array of SABR-calibrated volatilities
        /// </summary>
        public double[] SABRSkew()
        {
            double[] skew = new double[strike.Length];
            for (int i = 0; i < skew.Length; i++)
            {
                skew[i] = explicitAlphaVolatility(forward + shift, strike[i] + shift, expiry, alpha, beta, nu, rho);
            }
            return skew;
        }
        /// <summary>
        /// solve alpha coefficient
        /// </summary>
        private double solveForAlpha(double f, double t, double ATMVol, double beta, double nu, double rho)
        {
            double constant = -ATMVol * Math.Pow(f, (1 - beta));
            double linear = 1.0 + (2.0 - 3.0 * Math.Pow(rho, 2.0)) / 24.0 * Math.Pow(nu, 2.0) * t;
            double quadratic = 0.25 * rho * nu * beta * t / Math.Pow(f, (1.0 - beta));
            double cubic = Math.Pow(1.0 - beta, 2.0) * t / (24.0 * Math.Pow(f, (2.0 - 2.0 * beta)));
            //
            double root = 0.0;
            if (cubic > 0.0) root = MathTools.SmallestPositiveCubicRoot(cubic, quadratic, linear, constant);
            if (cubic == 0.0) root = MathTools.SmallestPositiveQuadraticRoot(quadratic, linear, constant);
            return root;
        }
        /// <summary>
        /// calculate volatility with explicitly given alpha
        /// </summary>
        private double explicitAlphaVolatility(double f, double x, double t, double alpha, double beta, double nu, double rho)
        {
            double[] SABR = new double[3];
            double sabrz, y;
            //
            // U-term
            SABR[0] = alpha / (Math.Pow(f * x, (1.0 - beta) / 2.0))
                * (1.0 + Math.Pow(1.0 - beta, 2.0) / 24.0 * Math.Pow(Math.Log(f / x), 2.0)
                + (Math.Pow(1.0 - beta, 4.0) / 1920.0) * Math.Pow(Math.Log(f / x), 4.0));
            //
            if (Math.Abs(f - x) > Math.Pow(10, -8.0))
            {
                sabrz = (nu / alpha) * Math.Pow(f * x, (1 - beta) / 2.0) * Math.Log(f / x);
                y = (Math.Sqrt(1.0 - 2.0 * rho * sabrz + Math.Pow(sabrz, 2.0)) + sabrz - rho) / (1.0 - rho);
                //
                if (Math.Abs(y - 1.0) < Math.Pow(10, -8.0))
                {
                    SABR[1] = 1.0;
                }
                else
                {
                    if (y > 0.0)
                    {
                        SABR[1] = sabrz / Math.Log(y);
                    }
                    else
                    {
                        SABR[1] = 1.0;
                    }
                }
            }
            else
            {
                SABR[1] = 1.0;
            }
            //
            // V-term
            SABR[2] = 1.0 + (((Math.Pow(1.0 - beta, 2.0) / 24.0) * (Math.Pow(alpha, 2.0) / (Math.Pow(f * x, 1.0 - beta))))
                + ((rho * beta * nu * alpha) / (4.0 * (Math.Pow(f * x, (1.0 - beta) / 2.0))))
                + ((((2.0 - 3.0 * Math.Pow(rho, 2.0)) / 24.0) * Math.Pow(nu, 2.0)))) * t;
            //
            return SABR[0] * SABR[1] * SABR[2];
        }
        /// <summary>
        /// calculate volatility without explicitly given alpha
        /// </summary>
        private double implicitAlphaVolatility(double f, double x, double t, double vol, double beta, double nu, double rho)
        {
            this.alpha = solveForAlpha(f, t, vol, beta, nu, rho);
            return explicitAlphaVolatility(f, x, t, alpha, beta, nu, rho);
        }
        /// <summary>
        /// callback method for LMA solver in order to calculate vector of function values
        /// </summary>
        void IFunctionVector.callback(double[] x, double[] fi, object obj)
        {
            // calculate squared differences of SABR volatility and market volatility 
            // assign calculated squared differences into function vector fi
            for (int i = 0; i < fi.Count(); i++)
            {
                double SABR_VOL = 0.0;
                //
                // solve for three coefficients : alpha = x[0], rho = x[1] and nu = x[2]
                if (x.Count() == 3) SABR_VOL = explicitAlphaVolatility(forward + shift, strike[i] + shift, expiry, x[0], beta, x[2], x[1]);
                //
                // solve for two coefficients (alpha is implicitly solved) : rho = x[0] and nu = x[1]
                if (x.Count() == 2) SABR_VOL = implicitAlphaVolatility(forward + shift, strike[i] + shift, expiry, atmVolatility, beta, x[1], x[0]);
                fi[i] = Math.Pow(marketVolatility[i] - SABR_VOL, 2.0);
            }
            //
            // update class member data coefficients
            if (x.Count() == 3)
            {
                this.alpha = x[0];
                this.rho = x[1];
                this.nu = x[2];
            }
            else
            {
                this.rho = x[0];
                this.nu = x[1];
                this.alpha = solveForAlpha(forward + shift, expiry, atmVolatility, beta, this.nu, this.rho);
            }
        }
    }
    //
    //
    //
    //
    public static class MathTools
    {
        /// <summary>
        /// find the smallest positive root for a given cubic polynomial function
        /// </summary>
        /// <param name="cubic">coefficient of cubic term</param>
        /// <param name="quadratic">coefficient of quadratic term</param>
        /// <param name="linear">coefficient of linear term</param>
        /// <param name="constant">constant term</param>
        public static double SmallestPositiveCubicRoot(double cubic, double quadratic, double linear, double constant)
        {
            double[] roots = new double[3];
            double a, b, c, r, q, capA, capB, theta;
            double root = 0.0;
            //
            a = quadratic / cubic;
            b = linear / cubic;
            c = constant / cubic;
            q = (Math.Pow(a, 2.0) - 3.0 * b) / 9.0;
            r = (2.0 * Math.Pow(a, 3.0) - 9.0 * a * b + 27.0 * c) / 54.0;
            //   
            if ((Math.Pow(r, 2.0) - Math.Pow(q, 3.0)) >= 0.0)
            {
                capA = -Math.Sign(r) * Math.Pow(Math.Abs(r) + Math.Sqrt(Math.Pow(r, 2.0) - Math.Pow(q, 3.0)), (1.0 / 3.0));
                if (capA == 0.0)
                {
                    capB = 0.0;
                }
                else
                {
                    capB = q / capA;
                }
                root = capA + capB - a / 3.0;
            }
            else
            {
                theta = Math.Acos(r / Math.Pow(q, 1.5));
                //
                // three available roots
                roots[0] = -2.0 * Math.Sqrt(q) * Math.Cos(theta / 3.0) - a / 3.0;
                roots[1] = -2.0 * Math.Sqrt(q) * Math.Cos(theta / 3.0 + 2.0943951023932) - a / 3.0;
                roots[2] = -2.0 * Math.Sqrt(q) * Math.Cos(theta / 3.0 - 2.0943951023932) - a / 3.0;
                //
                // identify the smallest positive root
                if (roots[0] > 0.0) root = roots[0];
                if (roots[1] > 0.0) root = roots[1];
                if (roots[2] > 0.0) root = roots[2];
                if (roots[1] > 0.0 && roots[1] < root) root = roots[1];
                if (roots[2] > 0.0 && roots[2] < root) root = roots[2];
            }
            return root;
        }
        //
        /// <summary>
        /// find the smallest positive root for a given quadratic polynomial function
        /// throws exception in the case of finding complex roots
        /// </summary>
        /// <param name="quadratic">coefficient of quadratic term</param>
        /// <param name="linear">coefficient of linear term</param>
        /// <param name="constant">constant term</param>
        public static double SmallestPositiveQuadraticRoot(double quadratic, double linear, double constant)
        {
            double x1, x2;
            double root = 0.0;
            double discriminant = Math.Pow(linear, 2.0) - 4.0 * quadratic * constant;
            //
            // case : complex roots (throw exception)
            if (discriminant < 0.0) throw new Exception("No real roots available");
            //
            // case : one real root
            if (discriminant == 0.0) root = -linear / (2.0 * quadratic);
            //
            // case : two real roots
            if (discriminant > 0.0)
            {
                x1 = (-linear + Math.Sqrt(discriminant)) / (2.0 * quadratic);
                x2 = (-linear - Math.Sqrt(discriminant)) / (2.0 * quadratic);
                root = x1;
                if ((x2 < x1) && (x2 > 0.0)) root = x2;
            }
            return root;
        }
    }
}

Finally, thanks for reading my blog.
-Mike