Re: Model calibration with external optimizers

Posted by jarkki on
URL: http://quantlib.414.s1.nabble.com/Model-calibration-with-external-optimizers-tp7186p7188.html

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;
        }
       

}

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