#include <ql/processes/hestonprocess.hpp>
#include <ql/models/equity/hestonmodel.hpp>
#include <ql/models/equity/hestonmodelhelper.hpp>
#include <ql/pricingengines/vanilla/analytichestonengine.hpp>
#include <ql/pricingengines/vanilla/mceuropeanhestonengine.hpp>
#include <ql/pricingengines/blackformula.hpp>
#include <ql/time/calendars/target.hpp>
#include <ql/time/calendars/nullcalendar.hpp>
#include <ql/time/daycounters/actual365fixed.hpp>
#include <ql/time/daycounters/actual360.hpp>
#include <ql/time/daycounters/actualactual.hpp>
#include <ql/termstructures/yieldcurves/zerocurve.hpp>
#include <ql/termstructures/yieldcurves/flatforward.hpp>
#include <ql/math/optimization/levenbergmarquardt.hpp>
#include <ql/time/period.hpp>

//todo remove this later
#include <ql/quantlib.hpp>

#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
#include <test-suite/utilities.hpp>
#include <iostream>
#include <iomanip>
#include <vector>
#include <fstream>


using namespace QuantLib;
using namespace std;

int main(int argc, char* argv[]) {
    /* this example is taken from A. Sepp
       Pricing European-Style Options under Jump Diffusion Processes
       with Stochstic Volatility: Applications of Fourier Transform
       http://math.ut.ee/~spartak/papers/stochjumpvols.pdf
    */

	cout << "Testing Heston model calibration using DAX volatility data..." << endl;

    Date settlementDate(5, July, 2002); // this is the settlement date for the 
										  // original example
	//Date settlementDate = Date::todaysDate();

    Settings::instance().evaluationDate() = settlementDate;

    DayCounter dayCounter = Actual365Fixed();
    Calendar calendar = TARGET();

    Integer t[] = { 13, 41, 75, 165, 256, 345, 524, 703 };
    Rate r[] = { 0.0357,0.0349,0.0341,0.0355,0.0359,0.0368,0.0386,0.0401 };

    std::vector<Date> dates;
    std::vector<Rate> rates;
	std::vector<Rate> dividends;
    dates.push_back(settlementDate);
    rates.push_back(0.0357);
    Size i;
    for (i = 0; i < 8; ++i) {
        dates.push_back(settlementDate + t[i]);
        rates.push_back(r[i]);
		dividends.push_back(0.1);
    }

    Handle<YieldTermStructure> riskFreeTS(
                       boost::shared_ptr<YieldTermStructure>(
                                    new ZeroCurve(dates, rates, dayCounter)));

    //Handle<YieldTermStructure> dividendTS(
    //                               flatRate(settlementDate, 0.0, dayCounter));

	Handle<YieldTermStructure> dividendTS(
  			    boost::shared_ptr<YieldTermStructure>(
  			    new ZeroCurve(dates, dividends, dayCounter))); 

    Volatility v[] =
      { 0.6625,0.4875,0.4204,0.3667,0.3431,0.3267,0.3121,0.3121,
        0.6007,0.4543,0.3967,0.3511,0.3279,0.3154,0.2984,0.2921,
        0.5084,0.4221,0.3718,0.3327,0.3155,0.3027,0.2919,0.2889,
        0.4541,0.3869,0.3492,0.3149,0.2963,0.2926,0.2819,0.2800,
        0.4060,0.3607,0.3330,0.2999,0.2887,0.2811,0.2751,0.2775,
        0.3726,0.3396,0.3108,0.2781,0.2788,0.2722,0.2661,0.2686,
        0.3550,0.3277,0.3012,0.2781,0.2781,0.2661,0.2661,0.2681,
        0.3428,0.3209,0.2958,0.2740,0.2688,0.2627,0.2580,0.2620,
        0.3302,0.3062,0.2799,0.2631,0.2573,0.2533,0.2504,0.2544,
        0.3343,0.2959,0.2705,0.2540,0.2504,0.2464,0.2448,0.2462,
        0.3460,0.2845,0.2624,0.2463,0.2425,0.2385,0.2373,0.2422,
        0.3857,0.2860,0.2578,0.2399,0.2357,0.2327,0.2312,0.2351,
        0.3976,0.2860,0.2607,0.2356,0.2297,0.2268,0.2241,0.2320 };

    Handle<Quote> s0(boost::shared_ptr<Quote>(new SimpleQuote(4468.17)));
    Real strike[] = { 3400,3600,3800,4000,4200,4400,
                      4500,4600,4800,5000,5200,5400,5600 };

    std::vector<boost::shared_ptr<CalibrationHelper> > options;

    for (Size s = 0; s < 13; ++s) {
        for (Size m = 0; m < 8; ++m) {
            Handle<Quote> vol(boost::shared_ptr<Quote>(
                                                  new SimpleQuote(v[s*8+m])));

            Period maturity((int)((t[m]+3)/7.), Weeks); // round to weeks
            options.push_back(boost::shared_ptr<CalibrationHelper>(
                        new HestonModelHelper(maturity, calendar,
                                              s0->value(), strike[s], vol,
                                              riskFreeTS, dividendTS, true)));
        }
    }

	//initial values for the calibration
    Real v0=0.1;
    Real kappa=1.0;
    Real theta=0.1;
    Real sigma=0.5;
    Real rho=-0.5;

    boost::shared_ptr<HestonProcess> hestonProcess(new HestonProcess(
              riskFreeTS, dividendTS, s0, v0, kappa, theta, sigma, rho));

    boost::shared_ptr<HestonModel> hestonModel(new HestonModel(hestonProcess));

    boost::shared_ptr<PricingEngine> hestonEngine(
							new AnalyticHestonEngine(hestonModel, 64));

    for (i = 0; i < options.size(); ++i)
        options[i]->setPricingEngine(hestonEngine);

    LevenbergMarquardt om(1e-8, 1e-8, 1e-8);
    hestonModel->calibrate(options, om, EndCriteria(400, 40, 1.0e-8, 1.0e-8, 1.0e-8));

    Real sse = 0;
    for (i = 0; i < 13*8; ++i) {
        const Real diff = options[i]->calibrationError()*100.0;
        sse += diff*diff;
    }

    Real expected = 177.2; //see article by A. Sepp.
    
	if (std::fabs(sse - expected) > 1.0) {
        cout << "Failed to reproduce calibration error"
                   << "\n    calculated: " << sse
                   << "\n    expected:   " << expected << endl;
    }

	//parameters of the calibrated model
	cout << "theta: " << hestonModel->theta() << endl;
	cout << "kappa: " << hestonModel->kappa() << endl;
	cout << "sigma: " << hestonModel->sigma() << endl;
	cout << "rho  : " << hestonModel->rho() << endl;
	cout << "v0   : " << hestonModel->v0() << endl;

	//Obtain the volatility surface using the calculated
	//heston parameters

	v0 = hestonModel->v0();	
    kappa = hestonModel->kappa();
    theta = hestonModel->theta();
    sigma = hestonModel->sigma();
    rho = hestonModel->rho();

	//to price using heston
	boost::shared_ptr<HestonProcess> hestonPricingProcess(new HestonProcess(
                   riskFreeTS, dividendTS, s0, v0, kappa, theta, sigma, rho));

	boost::shared_ptr<PricingEngine> hestonPricingEngine(new AnalyticHestonEngine(
              boost::shared_ptr<HestonModel>(new HestonModel(hestonPricingProcess)), 144));


	//to get the implied volatility using BS
	Volatility volatility = v[0];  //TODO CHECK

	Handle<BlackVolTermStructure> flatVolTS(
			    boost::shared_ptr<BlackVolTermStructure>(
				    new BlackConstantVol(settlementDate, volatility, dayCounter)));

	boost::shared_ptr<StochasticProcess> stochasticProcess(
                 new BlackScholesMertonProcess(s0, dividendTS, riskFreeTS, flatVolTS));
                                    

	boost::shared_ptr<PricingEngine> europeanPricingEngine = 
		boost::shared_ptr<PricingEngine>(new AnalyticEuropeanEngine);

	ofstream outfile("implvol.txt");

	for (Size s = 0; s < 13; ++s) {
		for (Size m = 0; m < 8; ++m) {	
			Period maturity(((t[m]+3)/7.), Weeks); 
			Date exerciseDate = settlementDate + maturity.length()*Weeks;

		    boost::shared_ptr<StrikedTypePayoff> payoff(
                                     new PlainVanillaPayoff(Option::Put, strike[s]));

			boost::shared_ptr<Exercise> exercise(new EuropeanExercise(exerciseDate));

			VanillaOption option4pricing(hestonPricingProcess, payoff, exercise);

			option4pricing.setPricingEngine(hestonPricingEngine);
			Real price = option4pricing.NPV();
			//cout << fixed << setprecision(2) << price << "\t";
	
			boost::shared_ptr<Exercise> europeanExercise(
                                         new EuropeanExercise(exerciseDate));

			

			boost::shared_ptr<VanillaOption> option4implvol(
                  new EuropeanOption(stochasticProcess, payoff, europeanExercise, europeanPricingEngine));
		
			Volatility implVol = 0.0; // just to remove a warning...
			if (price > 0.0) {
				Size maxEvaluations = 100;
				Real tolerance = 1.0e-6;

				try {
			
					implVol = option4implvol->impliedVolatility(price, tolerance, maxEvaluations);
					cout << fixed << setprecision(3) << implVol << "\t";
					outfile << fixed << setprecision(3) << implVol << "\t";
				} catch (std::exception& e) {	
						cout << e.what() << endl;
				}
			}
			else { 
				cout << fixed << setprecision(3) << 0.0 << "\t";
				outfile << fixed << setprecision(3) << 0.0 << "\t";
			}
		}
		cout << endl;
		outfile << endl;
	}
	outfile.close();
}