Monday, December 26, 2016

XLW : Interfacing C++11 PathGenerator to Excel

In my previous post on generating paths for different types of one-factor processes, I was writing all processed paths sequentially into separate CSV files. Later, these files could be opened in Excel for any further use. This approach is far from being an efficient way to interface C++ program to Excel. While being possible to create in-house tool for this task, the actual implementation part requires relatively deep understanding and hands-on experience on Excel/C API. The amount of issues to be learned in order to produce truly generic and reliable tool, is far from being something which could be internalized in very short period of time. There are limited amount of books available on this topic, but the most recommended book is the one written by Steve Dalton.

This time, I wanted to present kind of "industry standard" way to accomplish this task using "easy-to-use" Excel/C++ interfacing tool, which has been there already since 2002. XLW is an application which wraps Excel/C API into a simple C++ interface, which can then be used to customize Excel with user-defined worksheet functions. For any newcomer into this issue, it is highly recommended first to watch instructional "how-to-use" clips from project homepage. The complete workflow of using XLW wrapper is presented there from downloading to debugging. Also, AdaptiveRisk blog is presenting extremely useful stuff, well enough to get started.


cppinterface.h


For this project, I have extracted a new XLW xll template and opened corresponding solution in my Visual Studio Express 2013. The content of this header file in my current project is presented below. I have method declarations for two functions, which are both returning matrix object (XLW data structure). It is my intention to use my previous PathGenerator project, in order to fill these matrix objects with paths using desired one-factor processes.

#ifndef TEST_H
#define TEST_H
//
#include "xlw/MyContainers.h"
#include <xlw/CellMatrix.h>
#include <xlw/DoubleOrNothing.h>
#include <xlw/ArgList.h>
#include <xlw/XlfServices.h>
//
using namespace xlw;
//
//<xlw:libraryname=XLWPathGenerator
//
// method for requesting vasicek paths
MyMatrix // return matrix of random paths following Vasicek SDE
//<xlw:volatile
GetPaths_Vasicek(double t, // time to maturity
double r, // current short rate
double longTermRate, // long-term average rate
double meanReversion, // mean reversion speed
double rateVolatility // rate volatility
);
//
// method for requesting GBM paths
MyMatrix // return matrix of random paths following Geometric Brownian Motion SDE
//<xlw:volatile
GetPaths_BrownianMotion(double t, // time to maturity
double s, // current spot rate
double rate, // risk-free rate
double volatility // volatility
);
//
#endif

Commenting may seem a bit strange first, but the following screenshot containing Excel function argument input box may help to catch the point.















source.cpp


Implementations for two methods declared in header file are presented below. Information concerning number of time steps (for a path) and number of paths to be created are extracted from matrix dimensions using XlfServices object. After this, desired OneFactorProcess and PathGenerator objects are created. Finally, PathGenerator object is used to process a path, which will be imported into resulting matrix object (paths) and returned for the client (Excel).

#include <cppinterface.h>
#include "PathGenerator.h"
#pragma warning (disable : 4996)
//
MyMatrix GetPaths_Vasicek(double t, double r, double longTermRate, double meanReversion, double rateVolatility)
{
 // request dimensions for calling matrix
 const unsigned int nPaths = XlfServices.Information.GetCallingCell().columns();
 const unsigned int nSteps = XlfServices.Information.GetCallingCell().rows();
 // create container for all processed paths
 MyMatrix paths(nSteps, nPaths);
 // create container for a single path to be processed
 MyArray path(nSteps);
 //
 // create vasicek process and path generator
 std::shared_ptr<MJProcess::OneFactorProcess> vasicek =
  std::shared_ptr<MJProcess::Vasicek>(new MJProcess::Vasicek(meanReversion, longTermRate, rateVolatility));
 PathGenerator<> shortRateProcess(r, t, vasicek);
 //
 // process paths using path generator
 for (unsigned int i = 0; i != nPaths; ++i)
 {
  shortRateProcess(path);
  // import processed path into paths container
  for (unsigned j = 0; j != nSteps; ++j)
  {
   paths[j][i] = path[j];
  }
 }
 return paths;
}
//
MyMatrix GetPaths_BrownianMotion(double t, double s, double rate, double volatility)
{
 // request dimensions for calling matrix
 const unsigned int nPaths = XlfServices.Information.GetCallingCell().columns();
 const unsigned int nSteps = XlfServices.Information.GetCallingCell().rows();
 // create container for all processed paths
 MyMatrix paths(nSteps, nPaths);
 // create container for a single path to be processed
 MyArray path(nSteps);
 //
 // create geometric brownian motion process and path generator
 std::shared_ptr<MJProcess::OneFactorProcess> brownianMotion =
  std::shared_ptr<MJProcess::GBM>(new MJProcess::GBM(rate, volatility));
 PathGenerator<> equityPriceProcess(s, t, brownianMotion);
 //
 // process paths using path generator
 for (unsigned int i = 0; i != nPaths; ++i)
 {
  equityPriceProcess(path);
  // import processed path into paths container
  for (unsigned j = 0; j != nSteps; ++j)
  {
   paths[j][i] = path[j];
  }
 }
 return paths;
}
//


In order to get this thing up and running, header file for PathGenerator has to be included. I have set the current XLL project as startup project. As a side, I have opened my PathGenerator project (containing header files for RandomGenerator, OneFactorProcess and PathGenerator). Since this side project is still unaccessible, it has to be linked to my current XLL project : Project - Properties - Configuration Properties - C/C++ - General - Additional Include Directories (Browse folder containing source files for the project to be linked). After completing these steps and building this project succesfully, I am finally ready to test the provided functionality in Excel.


Excel


After opening a new Excel, I need to drag-and-drop (or open from Excel) newly created xll template from my \\Projects\XLWTester\Debug folder to Excel. In my Excel (screenshot below), I have two boxes for input parameters (Vasicek, Brownian Motion) and two ranges for resulting one-factor process paths (36 steps, 15 paths). As soon as I hit F9 button, my one-factor paths will be re-created. Finally, it should be noted that the functions are both array formulas.





















Finally, thanks for reading this blog. Pleasant waiting for a new year for everybody.
-Mike

Sunday, December 18, 2016

C++11 : PathGenerator Class

This blog posting is a sequel for the post I published just a couple of days ago. Having tools for generating random numbers is nice, but I still wanted to put that Random Generator template class to do something even more useful, like to produce paths for asset prices.

The idea was to create some practical means for transforming random number paths into asset price paths, following any desired (one-factor) stochastic process. For this purpose, I have created one template class (PathGenerator), which is technically just a wrapper for RandomGenerator template class and OneFactorProcess polymorphic class hierarchy. The purpose of RandomGenerator is to produce random numbers from desired probability distribution (usually that is standard normal) and the purpose of OneFactorProcess implementation is to provide information for PathGenerator on how to calculate drift and diffusion coefficients for a chosen stochastic process.


PATHS


For a marketing reasons, let us see the end product first. Simulated asset price paths for Geometric Brownian Motion and Vasicek processes are presented in Excel screenshot below. Test program (presented below) has been created in a way, that all processed asset price paths for a chosen stochastic process are exported into CSV file (which can then be imported into Excel for further investigation).

















ONE-FACTOR PROCESS


Abstract base class (OneFactorProcess) is technically just an interface, which provides practical means for a client for customizing drift and diffusion functions for different types of stochastic processes. I decided to implement polymorphic class hierarchy, since class is still pretty compact place for storing private member data and corresponding algorithms using that member data.

In the first draft, I was actually implementing drift and diffusion coefficient algorithms by using functions and lambdas, which (being initially created in main program) would then have been used inside PathGenerator object. It was technically working well, but from the viewpoint of possible end user (say, having member data and algorithm implementations in different files) it would have been quite a different story.

#pragma once
//
namespace MikeJuniperhillOneFactorProcessLibrary
{
 /// <summary>
 /// Abstract base class for all one-factor processes
 /// is technically just an interface, which provides practical means
 /// for customizing drift and diffusion functions for different types of processes.
 /// </summary>
 class OneFactorProcess
 {
 public:
  virtual double drift(double x, double t) = 0;
  virtual double diffusion(double x, double t) = 0;
 };
 //
 /// <summary>
 /// Implementation for Vasicek short-rate model.
 /// </summary>
 class Vasicek : public OneFactorProcess
 {
 public:
  Vasicek(double meanReversion, double longTermRate, double rateVolatility)
   : meanReversion(meanReversion), longTermRate(longTermRate), rateVolatility(rateVolatility) { }
  double drift(double x, double t) override { return meanReversion * (longTermRate - x); }
  double diffusion(double x, double t) override { return rateVolatility; }
 private:
  double meanReversion;
  double longTermRate;
  double rateVolatility;
 };
 //
 /// <summary>
 /// Implementation for Geometric Brownian Motion.
 /// </summary>
 class GBM : public OneFactorProcess
 {
 public:
  GBM(double rate, double volatility) : rate(rate), volatility(volatility) { }
  double drift(double x, double t) override { return rate * x; }
  double diffusion(double x, double t) override { return x * volatility; }
 private:
  double rate;
  double volatility;
 };
}


PATH GENERATOR


PathGenerator object is using RandomGenerator object for creating random numbers from standard normal probability distribution, by using default seeder for default uniform generator (Mersenne Twister). At some point, this class was having several different constructors for client-given seeder and client-given probability distribution. However, for the sake of clarity, I decided to remove all that optionality. In most of the cases, we want to simulate random numbers from standard normal distribution and Mersenne Twister generator still does the uniform part pretty well. It should be noted, that if such a customizing need sometimes arises, this class can be modified accordingly.

#pragma once
//
#include "RandomGenerator.h"
#include "OneFactorProcess.h"
namespace MJRandom = MikeJuniperhillRandomGeneratorTemplate;
namespace MJProcess = MikeJuniperhillOneFactorProcessLibrary;
//
template <typename Generator = std::mt19937, typename Distribution = std::normal_distribution<double>>
class PathGenerator
{
public:
 /// <summary>
 /// Constructor for PathGenerator template class.
 /// <para> - using default seeder from chrono library </para>
 /// <para> - using standard normal distribution </para>
 /// <para> - using mersenne twister uniform generator </para>
 /// </summary>
 PathGenerator(double spot, double maturity, std::shared_ptr<MJProcess::OneFactorProcess>& process)
  : spot(spot), maturity(maturity), process(process)
 {
  // create default seeder lambda function for random generator
   this->seeder = [](void) -> unsigned long { return static_cast<unsigned long>
    (std::chrono::steady_clock::now().time_since_epoch().count()); };
  // create random generator
  generator = std::unique_ptr<MJRandom::RandomGenerator<Generator, Distribution>>
   (new MJRandom::RandomGenerator<Generator, Distribution>(this->seeder));
 }
 /// <summary> 
 /// Functor, which fills auxiliary vector reference with asset prices following a given stochastic process.
 /// </summary>  
 void operator()(std::vector<double>& v)
 {
  // transform initialized vector into a path containing random numbers
  (*generator)(v);
  //
  double dt = maturity / (v.size() - 1);
  double dw = 0.0;
  double s = spot;
  double t = 0.0;
  v[0] = s; // 1st path element is always the current spot price
  //
  // transform random number vector into a path containing asset prices
  for (auto it = v.begin() + 1; it != v.end(); ++it)
  {
   t += dt;
   dw = (*it) * std::sqrt(dt);
   (*it) = s + (*process).drift(s, t) * dt + (*process).diffusion(s, t) * dw;
   s = (*it);
  }
 }
private:
 double spot;
 double maturity;
 std::function<unsigned long(void)> seeder;
 std::shared_ptr<MJProcess::OneFactorProcess> process;
 std::unique_ptr<MJRandom::RandomGenerator<Generator, Distribution>> generator;
};


TESTER


Presented test program is creating two different one-factor processes (Geometric Brownian Motion, Vasicek) and using PathGenerator for simulating asset price paths. All processed paths for the both cases will then be printed into CSV file for further investigations (Excel). I dare to say, that the process of generating asset price paths for one-factor process is easy and straightforward with PathGenerator class. For testing purposes, RandomGenerator header file should be included in the project.

#include "PathGenerator.h"
#include <fstream>
//
// printer : after a simulation iteration, append a path content into a file
// this inefficient idiom is not suitable for anything else but crude testing
void Print(const std::vector<double>& v, const std::string& filePathName)
{
 std::ofstream file(filePathName, std::ios_base::app);
 for (auto& element : v) { file << std::to_string(element) << ";"; }
 file << std::endl;
}
//
int main()
{
 // define and delete output CSV test files
 const std::string fileForVasicekProcess = "C:\\temp\\paths.vasicek.csv";
 const std::string fileForBrownianMotion = "C:\\temp\\paths.brownianMotion.csv";
 std::remove(fileForVasicekProcess.c_str());
 std::remove(fileForBrownianMotion.c_str());
 //
 // create auxiliary vector container for path processing
 double t = 3.0; // time to maturity
 int nPaths = 250; // number of simulated paths
 int nSteps = 1095; // number of time steps in one path (3 years * 365 days = 1095 steps)
 // the size of vector must be the number of desired time steps 
 // plus one, since the first vector item will be the current asset spot price
 std::vector<double> path(nSteps + 1);
 //
 // example 1 : create paths using vasicek process
 double r = 0.0095;
 double longTermRate = 0.05;
 double meanReversion = 0.2;
 double rateVolatility = 0.0075;
 std::shared_ptr<MJProcess::OneFactorProcess> vasicek =
  std::shared_ptr<MJProcess::Vasicek>(new MJProcess::Vasicek(meanReversion, longTermRate, rateVolatility));
 PathGenerator<> shortRateProcess(r, t, vasicek);
 for (int i = 0; i != nPaths; ++i)
 {
  shortRateProcess(path); // generate a path
  Print(path, fileForVasicekProcess); // print a path
 }
 //
 // example 2 : create paths using geometric brownian process
 double s = 100.0;
 double rate = 0.02;
 double volatility = 0.25;
 std::shared_ptr<MJProcess::OneFactorProcess> brownianMotion =
  std::shared_ptr<MJProcess::GBM>(new MJProcess::GBM(rate, volatility));
 PathGenerator<> equityPriceProcess(s, t, brownianMotion);
 for (int i = 0; i != nPaths; ++i)
 {
  equityPriceProcess(path);  // generate a path
  Print(path, fileForBrownianMotion); // print a path
 }
 return 0;
}

Finally, thanks for reading my blog and have a very pleasant waiting time for Christmas.
-Mike



Thursday, December 15, 2016

C++11 : Template Class for Random Generator

"Anyone who attempts to generate random numbers by deterministic means is, 
of course, living in a state of sin." - John Von Neumann

Despite of the fact that the quote above is still highly relevant, we mostly get the job done well enough using those pseudo-random generators for our Monte Carlo stuff. So, this time I wanted to present a wrapper for simulating (sinful) uniform random numbers, which are then mapped to a chosen probability distribution. The main goal was to construct a lightweight design, using C++11 uniform generators, probability distributions as well as some other extremely useful tools (function, lambda) for the task.


TEMPLATE DIET


Instead of implementing any class hierarchies, I tried to keep this design as lightweighted as possible by using template class having only two template data types to be defined by a client : Generator (uniform random generator) and Distribution (probability distribution for transforming uniformly distributed random numbers). In a nutshell, a client may use any random number engine available in C++11, for creating desired amount of uniform random numbers. Moreover, a client may use default probability distribution (Standard Normal) for transforming created uniform random numbers or additionally, use any desired probability distribution available in C++11. Available uniform random engines and probability distributions have been nicely covered in here.


ULTIMATE SIN


It should be a bit easier to sell something, if the final product is nicely presented. Random numbers processed by example program (presented later) for four different random generator implementations are presented in Excel screenshot below.




















HEADER FILE


Random generator template class is presented below. Note, that this implementation does not require any auxiliary libraries, since all the required stuff is already in C++11. For usability reasons, I decided to keep all the actual method implementations in header file (RandomGenerator.h).

#pragma once
#include <algorithm>
#include <chrono>
#include <cmath>
#include <functional>
#include <iostream>
#include <math.h>
#include <memory>
#include <random>
#include <string>
#include <vector>
//
namespace MikeJuniperhillRandomGeneratorTemplate
{
 template <typename Generator = std::mt19937, typename Distribution = std::normal_distribution<double>>
 /// <summary> 
 /// Template class for creating random number paths using
 /// Mersenne Twister as default uniform random generator and 
 /// Standard Normal (0.0, 1.0) as default probability distribution.
 /// </summary>  
 class RandomGenerator
 {
 public:
  /// <summary>
  /// Constructor for using default probability distribution 
  /// for transforming uniform random numbers.
  /// </summary>  
  RandomGenerator(const std::function<unsigned long(void)>& seeder)
   : seeder(seeder)
  {
   // construct lambda method for processing standard normal random number
   randomGenerator = [this](double x)-> double
   {
    x = distribution(uniformGenerator);
    return x;
   };
   // seed uniform random generator with a client-given algorithm
   uniformGenerator.seed(seeder());
  }
  /// <summary> 
  /// Constructor for using client-given probability distribution 
  /// for transforming uniform random numbers.
  /// </summary>  
  RandomGenerator(const std::function<unsigned long(void)>& seeder, const Distribution& distribution)
   // constructor delegation for initializing other required member data
   : RandomGenerator(seeder)
  {
   // assign client-given probability distribution
   this->distribution = distribution;
  }
  /// <summary> 
  /// Functor filling vector with random numbers mapped to chosen probability distribution.
  /// </summary>  
  void operator()(std::vector<double>& v)
  {
   std::transform(v.begin(), v.end(), v.begin(), randomGenerator);
  }
 private:
  std::function<unsigned long(void)> seeder;
  std::function<double(double)> randomGenerator;
  Generator uniformGenerator;
  Distribution distribution;
 };
}


TESTER FILE


Example program (tester.cpp) for testing random generator functionality is presented below. The process of generating (pseudo) random numbers from a probability distribution is easy and straightforward with this template wrapper class.

#include "RandomGenerator.h"
namespace MJ = MikeJuniperhillRandomGeneratorTemplate;
//
// printer for a random path
void Print(const std::string& message, const std::vector<double>& v)
{
 std::cout << message << std::endl;
 for (auto& element : v) std::cout << element << std::endl;
 std::cout << std::endl;
}
//
int main()
{
 // construct lambda method for seeding uniform random generator
 std::function<unsigned long(void)> seeder = 
  [](void) -> unsigned long { return static_cast<unsigned long>
  (std::chrono::steady_clock::now().time_since_epoch().count()); };
 //
 // create vector for a path to be filled with 20 random numbers
 // from desired probability distribution, using desired uniform random generator
 int nSteps = 20;
 std::vector<double> path(nSteps);
 //
 // path generator : mersenne twister, standard normal
 MJ::RandomGenerator<> gen_1(seeder);
 gen_1(path);
 Print("mt19937 : x~N(0.0, 1.0)", path);
 //
 // path generator : minst_rand, standard normal
 MJ::RandomGenerator<std::minstd_rand> gen_2(seeder);
 gen_2(path);
 Print("minstd_rand : x~N(0.0, 1.0)", path);
 //
 // path generator : mersenne twister, normal(112.5, 1984.0)
 std::normal_distribution<double> nd(112.5, 1984.0);
 MJ::RandomGenerator<> gen_3(seeder, nd);
 gen_3(path);
 Print("mt19937 : x~N(112.5, 1984.0)", path);
 //
 // path generator : minstd_rand, exponential(3.14)
 std::exponential_distribution<double> ed(3.14);
 MJ::RandomGenerator<std::minstd_rand, std::exponential_distribution<double>> gen_4(seeder, ed);
 gen_4(path);
 Print("minstd_rand : x~Exp(3.14)", path);
 //
 return 0;
}

Finally, thanks again for reading this blog.
-Mike

Thursday, October 27, 2016

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

Saturday, October 15, 2016

Alglib : Ho-Lee calibration in C#

A couple of years ago I published post on Ho-Lee short interest rate model calibration using Microsoft Solver Foundation (MSF). Solver is still available (somewhere) but not supported or developed any further, since Microsoft terminated that part of their product development (if I have understood the story correctly). As a personal note I have to say, that MSF was actually not the easiest package to use, from programming point of view.

Life goes on and very fortunately there are a lot of other alternatives available. This time I wanted to present Alglib, which is offering quite impressive package of numerical stuff already in their non-commercial (free of charge) edition. In optimization category, there is Levenberg-Marquardt algorithm (LMA) implementation available, which can be used for multivariate optimization tasks. Within this post, I am re-publishing short interest rate model calibration scheme, but only using Alglib LMA for performing the actual optimization routine.

In the case that Alglib is completely new package, it is recommended to check out some library documentation first, since it contains a lot of simple "Copy-Paste-Run"-type of examples, which will help to understand how the flow of information is processed. As a personal note I have to say, that time and effort needed to "getting it up and running" is extremely low. Just find Alglib download page, download proper C# version, create a new C# project and add reference to DLL file (alglibnet2.dll) included in download folder.

The outcome


Market prices of zero-coupon bonds, volatility (assumed to be estimated constant) and short rate are used, in order to solve time-dependent theta coefficients for Ho-Lee short interest rate model. When these coefficients have been solved, the model can be used to price different types of interest rate products. The original data and other necessary details have been already presented in here.

Alglib LMA package offers a lot of flexibility, concerning how the actual optimization task will be solved. Basically, there are three different schemes : V (using function vector only), VJ (using function vector and first order partial derivatives, known as Jacobian matrix) and FGH (using function vector, gradient and second order partial derivatives, known as Hessian matrix). Solved time-dependent theta coefficients for the first two schemes (V, VJ) are presented in the screenshot below.












Due to the usage of analytical first order partial derivatives instead of numerical equivalents (finite difference), Vector-Jacobian scheme (VJ) is a bit more accurate and computationally much faster. For calculating required partial derivatives for Jacobian/Hessian matrix, one might find this online tool quite handy.


The program


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

A few words on this program. First of all, Alglib LMA has been wrapped inside static LevenbergMarquardtSolver class, in order to provide kind of a generic way to use that solver. Basically, LMA requires callback method(s) for calculating values for vector of functions and/or Jacobian matrix. In this program, these callback methods are first defined as interfaces, which (the both in this program) are going to be implemented in HoLeeZeroCouponCalibration class. This class is also storing all the "source data" (maturities, zero-coupon bond prices, volatility and short rate) to be used, as LMA is processing its data (coefficients to be solved) through these callback methods. Note, that LevenbergMarquardtSolver class has absolutely no information about external world, except that it can receive callback methods defined in IFunctionVector and IJacobianMatrix interfaces. With this scheme, we can separate all financial algorithms and data from the actual optimization task.

using System;
using System.Linq;
//
// 'type definitions' for making life a bit more easier with long class names
using LM = AlgLibTester.LevenbergMarquardtSolver;
using HoLee = AlgLibTester.HoLeeZeroCouponCalibration;

namespace AlgLibTester
{
    class Program
    {
        static void Main(string[] args)
        {
            // short rate volatility, short rate, vector of maturities, vector of zero-coupon bond prices
            double sigma = 0.00039;
            double r0 = 0.00154;
            double[] t = new double[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 };
            double[] zero = new double[] { 0.9964, 0.9838, 0.9611, 0.9344, 0.9059, 0.8769, 0.8478, 0.8189, 0.7905, 0.7626 };
            //
            // create callback functions wrapped inside a class, which 
            // implements IFunctionVector and IJacobianMatrix interfaces
            HoLee callback = new HoLee(sigma, r0, t, zero);
            //
            // container for initial guesses for theta coefficients, which are going to be solved
            double[] theta = null;
            //
            // Example 1 :
            // use function vector only (using 'V' mode of the Levenberg-Marquardt optimizer)
            theta = new double[] { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
            LM.Solve(ref theta, 10, (IFunctionVector)callback);
            Console.WriteLine("Using function vectors only");
            Console.WriteLine("iterations : {0}", LM.report.iterationscount);
            theta.ToList<double>().ForEach(it => Console.WriteLine(it));
            Console.WriteLine("");
            //
            // Example 2 :
            // use function vector and jacobian matrix (using 'VJ' mode of the Levenberg-Marquardt optimizer)
            theta = new double[] { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; 
            LM.Solve(ref theta, 10, (IFunctionVector)callback, (IJacobianMatrix)callback);
            Console.WriteLine("Using function vectors and Jacobian matrix");
            Console.WriteLine("iterations : {0}", LM.report.iterationscount);
            theta.ToList<double>().ForEach(it => Console.WriteLine(it));
        }
    }
    //
    //
    //
    // static class, which is wrapper for Alglib Levenberg-Marquardt solver
    // in order to make the usage of this particular solver a bit more 'generic'
    public static class LevenbergMarquardtSolver
    {
        // alglib report is exposed to public
        public static alglib.minlmreport report;
        private static alglib.minlmstate state;
        //
        // create minimization solver using only function vector ('V' vector method)
        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);
        }
        // method overloading : create minimization solver using function vector and Jacobian matrix ('VJ' vector-jacobian method)
        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);
        }
    }
    //
    //
    //
    // interface : definition for function vector callback method, required by LevenbergMarquardtSolver class
    public interface IFunctionVector
    {
        void callback(double[] x, double[] fi, object obj);
    }
    //
    //
    //
    // interface : definition for jacobian matrix callback method, optionally required by LevenbergMarquardtSolver class
    public interface IJacobianMatrix
    {
        void callback(double[] x, double[] fi, double[,] jac, object obj);
    }
    //
    //
    //
    // Ho-Lee calibration class, which implements IFunctionVector and IJacobianMatrix interfaces
    public class HoLeeZeroCouponCalibration : IFunctionVector, IJacobianMatrix
    {
        // parameters, which are required in order to calculate
        // zero-coupon bond prices using Ho-Lee short rate model
        private double sigma;
        private double r0;
        private double[] t;
        private double[] zero;
        //
        public HoLeeZeroCouponCalibration(double sigma, double r0, double[] t, double[] zero)
        {
            this.sigma = sigma;
            this.r0 = r0;
            this.t = t;
            this.zero = zero;
        }
        // callback method for calculating vector of values for functions
        void IFunctionVector.callback(double[] x, double[] fi, object obj)
        {
            // calculate squared differences of Ho-Lee prices (using theta coefficients) 
            // and market prices and then assign these differences into function vector fi
            for (int i = 0; i < fi.Count(); i++)
            {
                double HOLEE_ZERO = Math.Exp(-Math.Pow(t[i], 2.0) * x[i] / 2.0
                    + Math.Pow(sigma, 2.0) * Math.Pow(t[i], 3.0) / 6.0 - r0 * t[i]);
                fi[i] = Math.Pow(zero[i] - HOLEE_ZERO, 2.0);
            }
        }
        // callback method for calculating partial derivatives for jacobian matrix
        void IJacobianMatrix.callback(double[] x, double[] fi, double[,] jac, object obj)
        {
            double HOLEE_ZERO = 0.0;
            // part 1 : calculate squared differences of Ho-Lee prices (using theta coefficients) 
            // and market prices, then assign these squared differences into function vector fi
            for (int i = 0; i < fi.Count(); i++)
            {
                HOLEE_ZERO = Math.Exp(-Math.Pow(t[i], 2.0) * x[i] / 2.0
                    + Math.Pow(sigma, 2.0) * Math.Pow(t[i], 3.0) / 6.0 - r0 * t[i]);
                fi[i] = Math.Pow(zero[i] - HOLEE_ZERO, 2.0);
            }           
            // part 2 : calculate all partial derivatives for Jacobian square matrix
            // loop through m functions
            for (int i = 0; i < fi.Count(); i++)
            {
                // loop through n theta coefficients
                for (int j = 0; j < x.Count(); j++)
                {
                    double derivative = 0.0;
                    // partial derivative is non-zero only for diagonal cases
                    if (i == j)
                    {
                        HOLEE_ZERO = Math.Exp(-Math.Pow(t[j], 2.0) * x[j] / 2.0
                            + Math.Pow(sigma, 2.0) * Math.Pow(t[j], 3.0) / 6.0 - r0 * t[j]);
                        derivative = Math.Pow(t[j], 2.0) * (zero[j] - HOLEE_ZERO) * HOLEE_ZERO;
                    }
                    jac[i, j] = derivative;
                }
            }
        }
    }
}

Finally, thanks for reading my blog.
-Mike