/*
 Copyright (C) 2003, 2004 Ferdinando Ametrano

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it under the
 terms of the QuantLib license.  You should have received a copy of the
 license along with this program; if not, please email quantlib-dev@lists.sf.net
 The license is also available online at http://quantlib.org/html/license.html

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#include "compound_option_pricing_engine.h"
#include<ql/quantlib.hpp>
#include <ql/PricingEngines/blackformula.hpp>

namespace QuantLib {

			class F
		{
			DiscountFactor dividendDiscount_;
			DiscountFactor riskFreeDiscount_;
			Real variance_;
			boost::shared_ptr<PlainVanillaPayoff> payoff_;
			Real compoundStrike_;

		public:
			F(
				DiscountFactor dividendDiscount,
				DiscountFactor riskFreeDiscount,
				Real variance,
				boost::shared_ptr<PlainVanillaPayoff> payoff,
				Real compoundStrike):dividendDiscount_(dividendDiscount),
									 riskFreeDiscount_(riskFreeDiscount),
									 variance_(variance),
									 payoff_(payoff),
									 compoundStrike_(compoundStrike){}
			Real operator ()(Real S) const {

				Real forward=S*dividendDiscount_/riskFreeDiscount_;

			 BlackFormula black(forward,riskFreeDiscount_,variance_,payoff_);

			 return (black.value()-compoundStrike_);
			}
		};

    void CompoundOptionPricingEngine::calculate() const {


        QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
                   "not an European Option");

		Option::Type=arguments_.payoff->type();

        Date exercise = arguments_.exercise->lastDate();

        boost::shared_ptr<PlainVanillaPayoff> payoff =
            boost::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);
        QL_REQUIRE(payoff, "non-plain payoff given");

        const boost::shared_ptr<BlackScholesProcess>& process =
            arguments_.blackScholesProcess;


        Volatility volatility =
            process->blackVolatility()->blackVol(exercise, payoff->strike());
        Real variance =
            process->blackVolatility()->blackVariance(exercise,
                                                      payoff->strike());
        DiscountFactor riskFreeDiscount =
            process->riskFreeRate()->discount(exercise);

        DayCounter rfdc  = process->riskFreeRate()->dayCounter();
        DayCounter divdc = process->dividendYield()->dayCounter();
        DayCounter voldc = process->blackVolatility()->dayCounter();

		 Time maturity = divdc.yearFraction(
            process->dividendYield()->referenceDate(), exercise);

		 Date compoundExpiration = arguments_.compoundExpiration;
		 Real compoundStrike=arguments_.compoundStrike;
		 Option::Type compoundType=arguments_.compoundType;

		 Time compoundMaturity = divdc.yearFraction(
			 process->dividendYield()->referenceDate(),compoundExpiration);

		 Real r=process->riskFreeRate()->zeroRate(exercise, rfdc,
                                              Continuous, NoFrequency);
		 Real q=process->dividendYield()->zeroRate(exercise, divdc,
                                               Continuous, NoFrequency);

		 Real corr=QL_SQRT(maturity/compoundMaturity);

		 Time deltaMaturity=compoundMaturity-maturity;

		 Real deltaRiskFreeDiscount=QL_EXP(-r*deltaMaturity);
		 Real deltaDividendDiscount = QL_EXP(-q*deltaMaturity);




	

		F f(deltaDividendDiscount,deltaRiskFreeDiscount,variance,payoff,compoundStrike);

		Bisection bisector;

		Real Sstar=bisector.solveImpl(f,.001);

		Real S=process->stateVariable()->value();

		Real X1=payoff->strike();

		Real a1=(QL_LOG(S/Sstar)+(r-q+0.5*variance)*maturity)/(QL_SQRT(variance*maturity));
		Real a2=a1-QL_SQRT(variance*maturity);

		Real b1=(QL_LOG(S/compoundMaturity)+(r-q+.5*variance)*compoundMaturity)/(QL_SQRT(variance*compoundMaturity));
		Real b2=b1-QL_SQRT(variance*compoundMaturity);

		BivariateCumulativeNormalDistribution Mcall(corr);
		BivariateCumulativeNormalDistribution Mput(-corr);


		if compoundType==Option::Call && 



   
    }
}