using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.Statistics;

namespace QuantConnect.Algorithm.CSharp
{
    /// <summary>
    /// This algorithm uses Math.NET Numerics library, specifically Linear Algebra object (Vector and Matrix) and operations, in order to solve a portfolio optimization problem.
    /// </summary>
    public class PortfolioOptimizationNumericsAlgorithm : QCAlgorithm
    {
        private string[] _symbols = new string[]
        {
        	// Using Meb Faber's GTAA paper assets:
            "SPY",   // 
    		"EFA",   // 
    		"TIP",   // 
    		"GSG",   // 
    		"VNQ"    //
    		// Find more symbols here: http://quantconnect.com/data
        };
        private const double _targetReturn = 0.1;
        private const double _riskFreeRate = 0.01;
        private double _lagrangeMultiplier;
        private double _portfolioRisk;
        private Matrix<double> Sigma;
        private List<SymbolData> SymbolDataList;

        public Vector<double> DiscountMeanVector
        {
            get
            {
                if (SymbolDataList == null)
                {
                    return null;
                }

                return 
                    Vector<double>.Build.DenseOfArray(SymbolDataList.Select(x => (double)x.Return).ToArray()) -
                    Vector<double>.Build.Dense(SymbolDataList.Count, _riskFreeRate);
            }
        }

        /// <summary>
        /// Initialise the data and resolution required, as well as the cash and start-end dates for your algorithm. All algorithms must initialized.
        /// </summary>
        public override void Initialize()
        {
            SetCash(100000);           //Set Strategy Cash
			SetStartDate(2013, 1, 1);  //Set Start Date
            SetEndDate(DateTime.Now.AddDays(-1));   //Set End Date
            
            SymbolDataList = new List<SymbolData>();

            foreach (var symbol in _symbols)
            {
                AddEquity(symbol, Resolution.Daily);
                SymbolDataList.Add(new SymbolData(symbol, History(symbol, 200, Resolution.Daily)));
            }

			Schedule.On(DateRules.MonthStart(), TimeRules.At(new TimeSpan(12, 0, 0)), () =>
            {
				ComputeWeights();

             	foreach (var symbolData in SymbolDataList.OrderBy(x => x.Weight))
                {
                    SetHoldings(symbolData.Symbol, symbolData.Weight);
                    Debug(Time.ToShortDateString() + " Purchased Stock: " + symbolData);
                }
            });
            
            //ComputePortfolioRisk();
        }

        /// <summary>
        /// OnData event is the primary entry point for your algorithm. Each new data point will be pumped in here.
        /// </summary>
        /// <param name="data">Slice object keyed by symbol containing the stock data</param>
        public override void OnData(Slice data)
        {
            foreach (var symbolData in SymbolDataList)
            {
            	if(data.ContainsKey(symbolData.Symbol))
            	{
                	symbolData.Update(data[symbolData.Symbol]);
            	}
            }
        }

        /// <summary>
        /// Computes Lagrange Multiplier
        /// </summary>
        private void ComputeLagrangeMultiplier()
        {
            var denominatorMatrix = DiscountMeanVector * Sigma.Inverse() * DiscountMeanVector.ToColumnMatrix();
			var denominator = denominatorMatrix.ToArray().First();
            _lagrangeMultiplier = denominator == 0 ? 0.0 : (_targetReturn - _riskFreeRate) / denominator;
        }

        /// <summary>
        /// Computes weight for each risky asset
        /// </summary>
        private void ComputeWeights()
        {
        	// Diagonal Matrix with each security risk (standard deviation)
            var S = Matrix<double>.Build.DenseOfDiagonalArray(SymbolDataList.Select(x => (double)x.Risk).ToArray());

            // Computes Correlation Matrix (using Math.NET Numerics Statistics)
            var allHistoryBars = new List<double[]>();
            SymbolDataList.ForEach(x => allHistoryBars.Add(x.History));
            var R = Correlation.PearsonMatrix(allHistoryBars);

            // Computes Covariance Matrix (using Math.NET Numerics Linear Algebra)
            Sigma = S * R * S;

            ComputeLagrangeMultiplier();
            
            var weights = _lagrangeMultiplier * Sigma.Inverse() * DiscountMeanVector.ToColumnMatrix();

            for (var i = 0; i < weights.RowCount; i++)
            {
                SymbolDataList[i].SetWeight(weights.ToArray()[i, 0]);
            }
        }

        /// <summary>
        /// Computes Portfolio Risk
        /// </summary>
        private void ComputePortfolioRisk()
        {
            var weights = Vector<double>.Build.DenseOfArray(SymbolDataList.Select(x => (double)x.Return).ToArray());
            var portfolioVarianceMatrix = weights * Sigma * weights.ToColumnMatrix();
            _portfolioRisk = Math.Sqrt(portfolioVarianceMatrix.ToArray().First());
            Log(string.Format("Lagrange Multiplier: {0,7:F4}", _lagrangeMultiplier));
            Log(string.Format("Portfolio Risk:      {0,7:P2} ", _portfolioRisk));
        }

        /// <summary>
        /// Symbol Data class to store security data (Return, Risk, Weight)
        /// </summary>
        class SymbolData
        {
            private RateOfChange _roc;
            private RollingWindow<double> _rollingHistory;
            private SimpleMovingAverage _sma;
            private StandardDeviation _std;
            public Symbol Symbol { get; private set; }
            public decimal Return { get { return _sma.Current; }  }
            public decimal Risk { get { return _std.Current; } }
            public decimal Weight { get; private set; }
            public double[] History
            {
                get
                {
                    return _rollingHistory.Select(x => x).ToArray();
                }
            }

            public SymbolData(Symbol symbol, IEnumerable<BaseData> history)
            {
                Symbol = symbol;
                Weight = 0m;
                _roc = new RateOfChange(2);
                _sma = new SimpleMovingAverage(200).Of(_roc);
                _std = new StandardDeviation(200).Of(_roc);
                _rollingHistory = new RollingWindow<double>(200);

                foreach (var data in history)
                {
                    Update(data);
                }
            }

            public void Update(BaseData data)
            {
            	if(data == null)
            	{
            		return;
            	}
            	else
            	{
                	_roc.Update(data.Time, data.Value);
                	_rollingHistory.Add((double)data.Value);
            	}
            }

            public void SetWeight(double value)
            {
            	Weight = value.IsNaNOrZero() ? 0m : (decimal)value;
            }

            public override string ToString()
            {
                return string.Format("{0}: {1,10:P2}\t{2,10:P2}\t{3,10:P2}", Symbol.Value, Weight, Return, Risk);
            }
        }
    }
}