namespace QuantConnect 
{   
    /*
    *   Example of how to use Put pricing with Black-Scholes in an algorithm
    *  		by: Jean-Paul van Brakel
    */
    public class BasicTemplateAlgorithm : QCAlgorithm
    {
    	// ticker to be used
    	private readonly string _ticker = "AAPL";
    	
    	// number of periods to be used in volatility calculation
    	private static int _vol_periods = 14;
    	private readonly RollingWindow<TradeBar> PriceHistory = new RollingWindow<TradeBar>(_vol_periods);
    	
    	// define option maturity date
    	private readonly DateTime _maturityDate = new DateTime(2015, 1, 16); // third friday of the month
    	
        //Initialize the data and resolution you require for your strategy:
        public override void Initialize() 
        {
            //Start and End Date range for the backtest:
            SetStartDate(2015, 1, 1); 
            SetWarmup(TimeSpan.FromDays(_vol_periods));
            SetEndDate(2015, 1, 2);
            
            //Cash allocation
            SetCash(25000);
            
            //Add as many securities as you like. All the data will be passed into the event handler:
            AddSecurity(SecurityType.Equity, _ticker, Resolution.Minute);
            
            //Initialise plot
			Chart plotter = new Chart("DerivativePlot", ChartType.Stacked);
	        plotter.AddSeries(new Series("Price", 		SeriesType.Line));
			plotter.AddSeries(new Series("Put price", 	SeriesType.Line));
			AddChart(plotter);
        }

        //Data Event Handler: New data arrives here.
        public void OnData(TradeBars data) 
        {   
        	PriceHistory.Add(data[_ticker]);
        	if (!PriceHistory.IsReady) return;
        	
        	// specify option settings here:
        	double price 		= (double)PriceHistory[0].Close;
        	double strike   	= 110;		// strike price of option
        	double rate 		= 0.05;		// risk-free rate of return to use in calculation
        	// recalculate annualised time to maturity
			double maturity 	= (BlackScholes.CountWeekDays(data.Time,_maturityDate)/250); 
			double[] _p_history = new double[PriceHistory.Count];
			for (int i = 0; i < PriceHistory.Count; i++)
				// copy close (you can change this to your liking)
				_p_history[i] = (double)PriceHistory[i].Close;
			// approximate volatility with historical volatility of the underlying
			double volatility 	= BlackScholes.HistoricalVolatility(_p_history); 
        	double yield 		= 0.01;		// approximation to the annualised dividend yields(%) for AAPL
        	
        	// calculate Black-Scholes option value for a European PUT (also approximation to American PUT)
        	double _optionPrice = BlackScholes.blsPut(price, strike, rate, maturity, volatility, yield);
        	
        	if (data.Time >= StartDate) {
				Plot("DerivativePlot", "Price", price);
            	Plot("DerivativePlot", "Put price", _optionPrice); 	
        	}
        	
        	// put your actual trading logic here:
        	if (!Portfolio.HoldStock) {
        		Order(_ticker, -100);
        	}
        }
    }
}

namespace QuantConnect {

    /// <summary>
    /// Black Scholes class to calculate the theoretical prices of European options (and Americans as an approximation)
	/// </summary>
    public static class BlackScholes
    {
    	private static double[] nsia = { 2.50662823884, -18.61500062529, 41.39119773534, -25.44106049637 };
        private static double[] nsib = { -8.4735109309, 23.08336743743, -21.06224101826, 3.13082909833 };
        private static double[] nsic = { 0.3374754822726147, 0.9761690190917186, 0.1607979714918209, 0.0276438810333863, 0.0038405729373609, 0.0003951896511919, 0.0000321767881768, 0.0000002888167364, 0.0000003960315187 };

        //cumulative normal distribution function
        private static double CND(double X)
        {
            double L = 0.0;
            double K = 0.0;
            double dCND = 0.0;
            const double a1 = 0.31938153;
            const double a2 = -0.356563782;
            const double a3 = 1.781477937;
            const double a4 = -1.821255978;
            const double a5 = 1.330274429;
            L = Math.Abs(X);
            K = 1.0 / (1.0 + 0.2316419 * L);
            dCND = 1.0 - 1.0 / Math.Sqrt(2 * Convert.ToDouble(Math.PI.ToString())) *
                Math.Exp(-L * L / 2.0) * (a1 * K + a2 * K * K + a3 * Math.Pow(K, 3.0) +
                a4 * Math.Pow(K, 4.0) + a5 * Math.Pow(K, 5.0));

            if (X < 0)
            {
                return 1.0 - dCND;
            }
            else
            {
                return dCND;
            }
        }

        //function phi
        private static double phi(double x)
        {
            double phi = 0.0;
            phi = Math.Exp(-x * x / 2) / Math.Sqrt(2 * Math.PI);
            return phi;
        }

        public static double NORMSINV(double probability)
        {
            double r = 0;
            double x = 0;
            x = probability - 0.5;
            if (Math.Abs(x) < 0.42)
            {
                r = x * x;
                r = x * (((nsia[3] * r + nsia[2]) * r + nsia[1]) * r + nsia[0]) / ((((nsib[3] * r + nsib[2]) * r + nsib[1]) * r + nsib[0]) * r + 1);
                return r;
            }
            r = probability;
            if (x > 0)
                r = 1 - probability;

            r = Math.Log(-Math.Log(r));
            r = nsic[0] + r * (nsic[1] + r * (nsic[2] + r * (nsic[3] + r * (nsic[4] + r * (nsic[5] + r * (nsic[6] + r * (nsic[7] + r * nsic[7])))))));
            return Math.Abs(r);
        }
        
        public static double NORMINV(double probability, double mean, double standard_deviation)
        {
            return (NORMSINV(probability) * standard_deviation + mean);
        }

        public static double NORMINV(double probability, double[] values)
        {
            return NORMINV(probability, Mean(values), StandardDeviation(values));
        }

        public static double Mean(double[] values)
        {
            double tot = 0;
            foreach (double val in values)
                tot += val;

            return (tot / values.Length);
        }
        
        public static double StandardDeviation(double[] values)
        {
            return Math.Sqrt(Variance(values));
        }
        
        public static double Variance(double[] values)
        {
            double m = Mean(values);
            double result = 0;
            foreach (double d in values)
                result += Math.Pow((d - m), 2);

            return (result / values.Length);
        }

		// calculate centered historical volatility over the last two weeks
		// derived from the annualized log difference of current price and last price
		public static double HistoricalVolatility(double[] historicalPrices) {
			double[] logDifferences = new double[historicalPrices.Length-1];
			for (int i = 1; i < historicalPrices.Length; i++) {
				logDifferences[i-1] = Math.Log(historicalPrices[i]/historicalPrices[i-1]);
			}
			double meanLogDifferences = Mean(logDifferences);
			double sumSquaredVariances = 0;
			for (int i = 0; i < logDifferences.Length; i++) {
				logDifferences[i] = Math.Pow(logDifferences[i] - meanLogDifferences,2);
				sumSquaredVariances += Math.Pow(logDifferences[i] - meanLogDifferences,2);
			}
			double sqrtAverageVariance = Math.Sqrt(sumSquaredVariances/logDifferences.Length);
			double annualizedVariance = sqrtAverageVariance*Math.Sqrt(252);
			
			return annualizedVariance;
		}


        //Call pricer
        public static double blsCall(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            double Call = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            Call = Price * Math.Exp(-Yield * Time) * CND(d1) - Strike * Math.Exp(-Rate * Time) * CND(d2);
            return Call;
        }

        //Put pricer
        public static double blsPut(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            double Put = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            Put = Strike * Math.Exp(-Rate * Time) * CND(-d2) - Price * Math.Exp(-Yield * Time) * CND(-d1);
            return Put;
        }

        //delta for Call
        public static double blsdeltaCall(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            return Math.Exp(-Yield * Time) * CND(d1);
        }

        //delta for Put
        public static double blsdeltaPut(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            return Math.Exp(-Yield * Time) * CND(d1) - 1;
        }

        //gamma is the same for Put and Call
        public static double blsgamma(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            return Math.Exp(-Yield * Time) * phi(d1) / (Price * Volatility * Math.Sqrt(Time));
        }

        //vega is the same for Put and Call
        public static double blsvega(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            return Price * Math.Exp(-Yield * Time) * phi(d1) * Math.Sqrt(Time);
        }

        //theta for Call
        public static double blsthetaCall(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            return -Math.Exp(-Yield * Time) *( Price * phi(d1) * Volatility / (2 * Math.Sqrt(Time))) - Rate * Strike * Math.Exp(-Rate * Time) * CND(d2) + Yield * Price * Math.Exp(-Yield * Time) * CND(d1);
        }

        //theta for Put
        public static double blsthetaPut(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            return -Math.Exp(-Yield * Time) *( Price * phi(d1) * Volatility / (2 * Math.Sqrt(Time))) + Rate * Strike * Math.Exp(-Rate * Time) * CND(-d2) - Yield * Price * Math.Exp(-Yield * Time) * CND(-d1);
        }

        //rho for Call
        public static double blsrhoCall(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            return Strike * Time * Math.Exp(-Rate * Time) * CND(d2);
        }

        //rho for Put
        public static double blsrhoPut(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            return -Strike * Time * Math.Exp(-Rate * Time) * CND(-d2);
        }

        //volga is the same for Call and Put
        public static double blsvolga(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            return Price * Math.Exp(-Yield * Time) * phi(d1) * Math.Sqrt(Time) * d1 * d2 / Volatility;

        }

        //vanna is the same for Call and Put
        public static double blsvanna(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            double vanna = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            vanna = -Math.Exp(-Yield * Time) * phi(d1) * d2 / Volatility;
            return vanna;
        }

        //charm for Call
        public static double blscharmCall(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            double charmC = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            charmC = -Yield * Math.Exp(-Yield * Time) * CND(d1) + Math.Exp(-Yield * Time) * phi(d1) * (2 * (Rate - Yield) * Time - d2 * Volatility * Math.Sqrt(Time)) / (2 * Time * Volatility * Math.Sqrt(Time));
            return charmC;
        }

        //charm for Put
        public static double blscharmPut(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            double charmP = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            charmP = Yield * Math.Exp(-Yield * Time) * CND(-d1) - Math.Exp(-Yield * Time) * phi(d1) * (2 * (Rate - Yield) * Time - d2 * Volatility * Math.Sqrt(Time)) / (2 * Time * Volatility * Math.Sqrt(Time));
            return charmP;
        }

        //color is the same for Call and Put
        public static double blscolor(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            double color = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            color = -Math.Exp(-Yield * Time) * (phi(d1) / (2 * Price * Time * Volatility * Math.Sqrt(Time))) * (2 * Yield * Time + 1 + (2 * (Rate - Yield) * Time - d2 * Volatility * Math.Sqrt(Time)) * d1 / (2 * Time * Volatility * Math.Sqrt(Time)));
            return color;
        }

        //dual delta for Call
        public static double blsdualdeltaCall(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            double ddelta = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            ddelta = -Math.Exp(-Rate * Time) * CND(d2);
            return ddelta;
        }

        //dual delta for Put
        public static double blsdualdeltaPut(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            double ddelta = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            ddelta = Math.Exp(-Rate * Time) * CND(-d2);
            return ddelta;
        }

        //dual gamma is the same for Call and Put
        public static double blsdualgamma(double Price, double Strike, double Rate, double Time, double Volatility, double Yield)
        {
            double d1 = 0.0;
            double d2 = 0.0;
            double dgamma = 0.0;
            d1 = (Math.Log(Price / Strike) + (Rate - Yield + Volatility * Volatility / 2.0) * Time) / (Volatility * Math.Sqrt(Time));
            d2 = d1 - Volatility * Math.Sqrt(Time);
            dgamma = Math.Exp(-Rate * Time) * phi(d2) / (Strike * Volatility * Math.Sqrt(Time));
            return dgamma;
        }
        
		/// <summary>
        /// Counts the number of week days (business days). Not so robust. Better to use real calendar.
		/// </summary>
        public static double CountWeekDays(DateTime d0, DateTime d1)
        {
            int ndays = 1 + Convert.ToInt32((d1 - d0).TotalDays);
            int nsaturdays = (ndays + Convert.ToInt32(d0.DayOfWeek)) / 7;
            return (double) (ndays - 2 * nsaturdays
                   - (d0.DayOfWeek == DayOfWeek.Sunday ? 1 : 0)
                   + (d1.DayOfWeek == DayOfWeek.Saturday ? 1 : 0));
        }
    }

}