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


Saturday, October 8, 2016

C# : managing global configurations with XML serializer

Managing configurations within any non-trivial program can easily turn to be a tricky issue. One might have a lot of different types of configuration parameters, which might need to be used as input parameters for several different classes. Moreover, there might be a need for using different set of configurations and one should be able to switch to another set as easily as possible.

This time I wanted to present a scheme, what I have been using for hosting global configuration parameters in some of my applications. Presented scheme solves three following issues :
  1. Configuration set used by program executable can easily be changed to another set
  2. There is no need to touch App.Config file
  3. All configuration properties are global, static and safe (read-only)
First of all, a set of configurations will be processed (read) from XML file using XML serializer. Assume I have following two configuration sets as XML files, one for testing and another for production :



















When running program executable for the both configuration files (path to folder in which the actual configuration XML file exists, will be given as a parameter for executable), we get the following console printouts :










When the program starts, path to configuration folder is captured (command line argument). Calling any class property of static Configurations class will automatically trigger static class constructor. Static class constructor will then use XML serializer (ObjectXMLHandler<T> class) for creating Configuration object instance (nested class inside static Configurations class). Finally, properties of newly created Configuration object will be assigned to properties of static Configurations class.

For testing this scheme, just replicate XML configuration files, create a new console application (GlobalConfigurations) and copyPaste the following code into CS-file.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
using System.Xml.Serialization;

namespace GlobalConfigurations
{
    public class Program
    {
        public static string configurationFilePathName;
        //
        static void Main(string[] args)
        {
            configurationFilePathName = args[0];
            //
            // print configurations
            ConfigurationTester configurationTester = new ConfigurationTester();
            configurationTester.Print();
        }
    }
    //
    //
    //
    /// <summary>
    /// Class for testing global static read-only configurations 
    /// </summary>
    public class ConfigurationTester
    {
        public void Print()
        {
            // print configurations
            Console.WriteLine("{0}, {1}, {2}, {3}", 
                Configurations.ConfigurationString, 
                Configurations.ConfigurationBoolean, 
                Configurations.ConfigurationLong, 
                Configurations.ConfigurationDouble);
        }
    }
    //
    //
    //
    /// <summary>
    /// Static class for hosting global configurations 
    /// </summary>
    public static class Configurations
    {
        // static read-only properties
        public static readonly string ConfigurationString;
        public static readonly bool ConfigurationBoolean;
        public static readonly long ConfigurationLong;
        public static readonly double ConfigurationDouble;
        //
        /// <summary>
        /// Static constructor will be called immediately 
        /// when static class property will be called. 
        /// </summary>
        static Configurations()
        {
            // create instance of ObjectXMLHandler for Configuration type
            ObjectXMLHandler<Configuration> handler = new ObjectXMLHandler<Configuration>();
            //
            // create Configuration objects from XML files
            string filePathName = Program.configurationFilePathName;
            List<Configuration> configuration = handler.Deserialize(filePathName);
            //
            // assign properties from inner configuration class into static read-only properties
            ConfigurationString = configuration[0].configurationString;
            ConfigurationBoolean = configuration[0].configurationBoolean;
            ConfigurationLong = configuration[0].configurationLong;
            ConfigurationDouble = configuration[0].configurationDouble;
        }
        /// <summary>
        /// Nested class for inner configurations. 
        /// </summary>
        public class Configuration
        {
            public string configurationString;
            public bool configurationBoolean;
            public long configurationLong;
            public double configurationDouble;
            //
            /// <summary>
            /// Default constructor is needed to create configuration parameters 
            /// </summary>
            public Configuration()
            {
                // no implementation
            }
        }
    }
    //
    //
    //
    /// <summary>
    /// Generic template class for handling conversion 
    /// between object and XML presentation.
    /// </summary>
    public class ObjectXMLHandler<T>
    {
        private XmlSerializer serializer = null;
        private Stream stream = null;
        //
        /// <summary>
        /// Parameter constructor for creating an instance. 
        /// </summary>
        public ObjectXMLHandler()
        {
            serializer = new XmlSerializer(typeof(T));
        }
        /// <summary>
        /// Convert a list of objects of type T into XML files. 
        /// </summary>
        public void Serialize(List<T> objects, List<string> fileNames, string folderPathName)
        {
            if (objects.Count != fileNames.Count)
                throw new Exception("objects.Count != fileNames.Count");
            //
            int counter = 0;
            foreach (T t in objects)
            {
                stream = File.Create(Path.Combine(folderPathName, String.Concat(fileNames[counter])));
                serializer.Serialize(stream, t);
                stream.Close();
                counter++;
            }
        }
        /// <summary>
        /// Convert XML files in specific folder into a list of objects of type T.
        /// </summary>
        public List<T> Deserialize(string folderPathName)
        {
            List<T> objects = new List<T>();
            foreach (string t in Directory.GetFiles(folderPathName))
            {
                stream = File.OpenRead(t);
                objects.Add((T)serializer.Deserialize(stream));
                stream.Close();
            }
            return objects;
        }
    }
}

Thanks a lot for using your precious time and reading my blog.
-Mike