Model calibration with external optimizers

[jarkki]

Hi,

I am trying to use the Differential Evolution algorithm and the Levenberg-Marquardt implementation of Manolis Lourakis to calibrate for example the Heston, Bates and Merton's jump diffusion models. What is the 'right' or the fastest way to implement the cost function?

I think it could work this way:
1. Create a vector of pointers to options and set the engine/model/process for each option
2. Pass this vector and other relevant data to the optimization routine/cost function
3. In the cost function, update the model parameters with the ones given by the optimization algorithm.
4. In the cost function, Call the option->NPV() for each of the options

I just don't know how to update the model parameters in the cost function. I suppose I should call the model->setParams(newParams), but how do I do this from the option object? Or should I use the CalibrationHelper object?

Thanks,
Jarkki

[Luigi Ballabio]

On Mon, 2009-02-09 at 02:02 -0800, jarkki wrote:
> I am trying to use the Differential Evolution algorithm and the
> Levenberg-Marquardt implementation of Manolis Lourakis to calibrate for
> example the Heston, Bates and Merton's jump diffusion models. What is the
> 'right' or the fastest way to implement the cost function?

You don't have to.  If your model inherits from CalibratedModel, the
cost function is already written for you (along the lines you described;
you were on the right track.)  For instance, you can already use the
HestonModel or BatesModel classes this way.

What you have to do is:

- on the optimization side: write your algorithm as a class derived from
OptimizationMethod and implement its minimize() method.  The method is
passed an instance of the Problem class, from which it takes the cost
function, a constraint, and the initial value.  With those pieces, it
can perform its job in any way you see fit; when it's done, it should
store the result into the passed Problem (by calling the setCurrentValue
and setFunctionValue methods) and return the reason why it finished
(success, failure to converge, errors...)  For an example, you can look
at the available Simplex or Levenberg-Marquardt implementation.

- on the model side: if you want to implement a model (other than the
available Heston and Bates) inherit it from CalibratedModel. The base
class stores the parameters to be calibrated and provides the plumbing.
In your model class, you'll have to implement the generateArguments()
method so that it uses the current value of the parameters to set up any
object required by the pricing engines (for instance, in the HestonModel
class, it builds a Heston process.)  You'll also have to write at least
a pricing engine using your model, and a calibration helper using the
engine; see AnalyticHestonEngine and HestonModelHelper for examples.

Once you have the two pieces, the code to tie them together is already
there; you can just call

model.calibrate(helpers, optimizer, ...)

which will create the cost function, feed it to your optimizer, and
calibrate the model.  You can see how calibration is set up in the
Heston-model test cases.

Let me know if anything is still not clear.

Luigi


--

Present to inform, not to impress; if you inform, you will impress.
-- Fred Brooks



------------------------------------------------------------------------------
Create and Deploy Rich Internet Apps outside the browser with Adobe(R)AIR(TM)
software. With Adobe AIR, Ajax developers can use existing skills and code to
build responsive, highly engaging applications that combine the power of local
resources and data with the reach of the web. Download the Adobe AIR SDK and
Ajax docs to start building applications today-http://p.sf.net/sfu/adobe-com
_______________________________________________
QuantLib-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-users

[jarkki]

Thank you very much for your answer.

I have the optimization routine set up as you described(I guess), but my optimizer fails to converge.

I am trying to calibrate the Heston model for just one maturity, and to the prices, not to the volatilities. I think the problem might lie in the set up of the calibration helper and term structure objects. I also tried with some reference data/prices, and in the cost function of my optimizer, the "problem->values(parameters)" gave wrong prices. I would appreciate, if you have time, if you could check my code for errors.

Just for testing, I tried the calibration with creating the termstrct/process/model/engine objects each evaluation time in the cost function, then calling the NPV(), and the algorithm converged. (obviously MEGA slow)

Thanks again,
Jarkki


Here is the code:

***************************************
main function:

        //Take data
        vector<double> X = strikes();
        vector<double> calls = calls();
        Volatility v[] = {0.3984, 0.3808, 0.3455, 0.3194, 0.3039, 0.2785, 0.2646, 0.2373, 0.2260, 0.2129, 0.2049};

        //Set variables
        Real spot = 6229.0;
        Rate r = 0.059;
        double T = 0.0767; // 28/365
        double divid = 0.0;

        //Initial values
        Real v0 = 0.01;
        Real kappa = 2.0;
        Real theta = 0.01;
        Real sigma = 0.225;
        Real rho = 0.0;
               
        // Set dates
        Calendar calendar = TARGET();
        Date evaluationDate(18,Feb, 2000);
        Date expiryDate(17, Mar, 2000);
        DayCounter dayCounter = Actual365Fixed();
        std::vector<Date> dates;
        std::vector<Rate> rates;

        for(int i = 0; i < X.size(); ++i)
        {
                dates.push_back(expiryDate);
                rates.push_back(r);
        }
       
        // The date at which pricing is to be performed.
        Settings::instance().evaluationDate() = evaluationDate;

        //Spot price handle
        Handle<Quote> s0(boost::shared_ptr<Quote>(new SimpleQuote(spot)));

        //Risk free rate
        boost::shared_ptr<Quote> riskFreeRate(new SimpleQuote(r));
    boost::shared_ptr<YieldTermStructure> flatRate(new FlatForward(0, NullCalendar(), Handle<Quote>(riskFreeRate), dayCounter));
    Handle<YieldTermStructure> riskFreeTS(flatRate);
       
        //Dividend yield (no dividends)
    Handle<YieldTermStructure> dividendTS(boost::shared_ptr<YieldTermStructure>(new FlatForward(evaluationDate, 0.0, dayCounter)));

        //Create Heston process
        boost::shared_ptr<HestonProcess> process(new HestonProcess(riskFreeTS, dividendTS, s0, v0, kappa, theta, sigma, rho));
       
        //Create Heston model
        boost::shared_ptr<HestonModel> model(new HestonModel(process));
   
        //Create an analytical pricing engine
        Size integrationOrder = 192;
    boost::shared_ptr<PricingEngine> engine(new AnalyticHestonEngine(model, integrationOrder));

        //Days to maturity
        Period maturity(dayCounter.dayCount(evaluationDate, expiryDate), Days);

        //Vector for calibration helpers
    vector<boost::shared_ptr<CalibrationHelper> > options;

        //Loop through strikes and create helper objects
        for(int i = 0; i < X.size(); ++i)
        {
                //Create quote for volatility
                Handle<Quote> vol(boost::shared_ptr<Quote>(new SimpleQuote(v[i])));
               
                //Create helper for Heston calibration
                boost::shared_ptr<CalibrationHelper> option(new HestonModelHelper(maturity, calendar, s0->value(), X[i], vol, riskFreeTS, dividendTS, false));
               
                //Set the pricing engine for the helper
                option->setPricingEngine(engine);

                options.push_back(option);
       
        }
       


        //Structure to hold data for the optimization procedure
        CalibrationData hestonCalibData;
        hestonCalibData.XCall = X;
        hestonCalibData.calls = calls;
        hestonCalibData.volas = volas;
       
        //Arguments for the Differential Evolution optimizer
        const int hestonNDim = 5;
        const int hestonPopSize = 100;
        const int hestonMaxGener = 2000;
        double hestonMin[5] = {0.001, 0.001, 0.001, 0.001, -1.0};
        double hestonMax[5] = {2.0, 10.0, 2.0, 2.5, 1.0};

        //Create the DE optimizer
        QLHestonSolver hestonSolver(hestonNDim,hestonPopSize, hestonCalibData);

        //Setup the solver
        hestonSolver.Setup(hestonMin, hestonMax, stBest1Exp, 0.9, 1.0, time(0));
               
        cout << "Started calibration..." << endl;
       
        //Calibrate
        model->calibrate(options, hestonSolver, EndCriteria(hestonMaxGener,40, 1.0e-8, 1.0e-8, 1.0e-8));

**************************************************************
Then here is the minimize function:
EndCriteria::Type QLHestonSolver::minimize(Problem & P, const EndCriteria & endCriteria)
{
        P.reset();
        this->problem = &P;
        this->Solve(endCriteria.maxIterations());
        double * res = this->Solution();
       
        QuantLib::Array results(5);
        for(int i = 0; i < 5; ++i)
        {
                results[i] = res[i];
        }
       
       
        EndCriteria::Type ecType = EndCriteria::None;
       
        P.setCurrentValue(results);
       
        return ecType;
}

****************************************
And the cost function of the DE:
double QLHestonSolver::EnergyFunction(double beta[], bool & bAtSolution)
{
        QuantLib::Array p(5);
        p[0] = beta[0];
        p[1] = beta[1];
        p[2] = beta[2];
        p[3] = beta[3];
        p[4] = beta[4];
                                       
        if(this->problem->constraint().test(p))
        {
                const QuantLib::Array & theorPrices = this->problem->values(p);

                double sse = 0.0;
                for(int i = 0; i < theorPrices.size(); ++i)
                {
                        double distance = this->data.calls[i] - theorPrices[i];
                        sse += (distance*distance);

                }

               return sse;

        }else
        {
                return 100000000000000000000.0;
        }
       

}

***********************************************************

[Luigi Ballabio]

On Thu, 2009-02-19 at 04:28 -0800, jarkki wrote:
> I have the optimization routine set up as you described(I guess), but my
> optimizer fails to converge.
>
> I am trying to calibrate the Heston model for just one maturity, and to the
> prices, not to the volatilities. I think the problem might lie in the set up
> of the calibration helper and term structure objects. I also tried with some
> reference data/prices, and in the cost function of my optimizer, the
> "problem->values(parameters)" gave wrong prices. I would appreciate, if you
> have time, if you could check my code for errors.

I just had a quick look at the setup (I didn't run the code) and there's
a problem.  (Incidentally, that's one that returns over and over and
keeps fooling people. We'll have to do something about it some time.)
Anyway, if you do:

> //Days to maturity
> Period maturity(dayCounter.dayCount(evaluationDate, expiryDate), Days);

and

> //Create helper for Heston calibration
> boost::shared_ptr<CalibrationHelper> option(new
> HestonModelHelper(maturity, calendar, s0->value(), X[i], vol, riskFreeTS,
> dividendTS, false));

You'll end up with the wrong expiration dates for the helpers, because
maturity (which is a period of n calendar days) will be interpreted by
the helper as a number of business days instead.  May you try replacing
the first by

Period maturity(calendar.businessDaysBetween(evaluationDate,
expiryDate), Days);

and see what happens?

Luigi


--

Green's Law of Debate:
Anything is possible if you don't know what you're talking about.



------------------------------------------------------------------------------
Open Source Business Conference (OSBC), March 24-25, 2009, San Francisco, CA
-OSBC tackles the biggest issue in open source: Open Sourcing the Enterprise
-Strategies to boost innovation and cut costs with open source participation
-Receive a $600 discount off the registration fee with the source code: SFAD
http://p.sf.net/sfu/XcvMzF8H
_______________________________________________
QuantLib-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-users