/*
 * MonteCarloSimulation.cpp
 *
 *  Created on: 20.05.2011
 *      Author: sven
 */

#include "Normaldistribution.h"
#include "MonteCarloSimulation.h"
#include <iomanip>
#include <iostream>
#include <algorithm>
#include <vector>
#include <fstream>

static const unsigned NUMBER_OF_SIMULATIONS = 10000;

const double PORTFOLIO[] = {200000.0,150000.0,250000.0,300000.0,200000.0};

MonteCarloSimulation::MonteCarloSimulation()
{
	for (unsigned i=0; i<NUMBER_OF_SIMULATIONS; i++) {
		mOneDaySimulatedCourseChanges.push_back(std::vector<double>(nenEndOfCompanies));
		mTenDaySimulatedCourseChanges.push_back(std::vector<double>(nenEndOfCompanies));
	}
}

MonteCarloSimulation::~MonteCarloSimulation() {
}

void MonteCarloSimulation::doSimulations(const Parametrisation& param, const CholeskyDecomposition& cholDecomp)
{
	std::vector<double> y(nenEndOfCompanies);
	std::vector<double> mu(nenEndOfCompanies);
	Normaldistribution norm;

	for (int i=0; i<2; i++) {
		enHoldingTime d = (i==0)?nenOneDay:nenTenDays;
		const CholeskyMatrix& D = cholDecomp.getMatrix(d);
		param.getArithAverages(mu,d);
		for (unsigned i=0; i< NUMBER_OF_SIMULATIONS; i++) {
			norm.calcRandomVector(D,mu,y);
			if (d == nenOneDay) {
				for (unsigned j=0; j<nenEndOfCompanies;j++) {
					mOneDaySimulatedCourseChanges.at(i).at(j) = y.at(j);
				}
			} else if (d == nenTenDays) {
				for (unsigned j=0; j<nenEndOfCompanies;j++) {
					mTenDaySimulatedCourseChanges.at(i).at(j) = y.at(j);
				}
			}
		}
	}
	std::vector<double> portfolio(nenEndOfCompanies);
	double portfolioValue = 0.0;
	for (unsigned c=0; c < nenEndOfCompanies; c++)
	{
		portfolio.at(c) = param.getLastValue(static_cast<enCompany>(c) ) * PORTFOLIO[c];
		portfolioValue += portfolio.at(c);
	}

	std::vector< std::vector<double> > simulatedPortfolioOneDay;
	std::vector< std::vector<double> > simulatedPortfolioTenDays;
	std::vector<double> simulatedPortfolioValuesOneDay(NUMBER_OF_SIMULATIONS);
	std::vector<double> simulatedPortfolioValuesTenDays(NUMBER_OF_SIMULATIONS);

	for (unsigned i=0; i<NUMBER_OF_SIMULATIONS; i++) {
		simulatedPortfolioOneDay.push_back(std::vector<double>(nenEndOfCompanies));
		simulatedPortfolioTenDays.push_back(std::vector<double>(nenEndOfCompanies));
		simulatedPortfolioValuesOneDay.at(i) = 0.0;
		simulatedPortfolioValuesTenDays.at(i) = 0.0;
	}

	for (unsigned i=0; i<NUMBER_OF_SIMULATIONS; i++) {
		for (unsigned j=0; j<nenEndOfCompanies;j++) {
			simulatedPortfolioOneDay.at(i).at(j) = portfolio.at(j) + (portfolio.at(j) * mOneDaySimulatedCourseChanges.at(i).at(j) / 100.0);
			simulatedPortfolioTenDays.at(i).at(j) = portfolio.at(j) + (portfolio.at(j) * mTenDaySimulatedCourseChanges.at(i).at(j) / 100.0);
			simulatedPortfolioValuesOneDay.at(i) += simulatedPortfolioOneDay.at(i).at(j);
			simulatedPortfolioValuesTenDays.at(i) += simulatedPortfolioTenDays.at(i).at(j);
		}
//		std::cout << std::setw(20) << std::fixed << (portfolioValue - simulatedPortfolioValuesOneDay.at(i))
//				<< std::setw(20) << (portfolioValue - simulatedPortfolioValuesTenDays.at(i))
//				<< std::endl;
	}
	std::sort(simulatedPortfolioValuesOneDay.begin(), simulatedPortfolioValuesOneDay.end());
	std::sort(simulatedPortfolioValuesTenDays.begin(), simulatedPortfolioValuesTenDays.end());

	std::vector<double> portfolioValueChangesOneDay(NUMBER_OF_SIMULATIONS);
	std::vector<double> portfolioValueChangesTenDays(NUMBER_OF_SIMULATIONS);
	std::ofstream of;
	of.open("VaR.csv",std::ios::out);
	for (unsigned i=0; i<NUMBER_OF_SIMULATIONS; i++) {
		portfolioValueChangesOneDay.at(i) = simulatedPortfolioValuesOneDay.at(i) - portfolioValue;
		portfolioValueChangesTenDays.at(i) = simulatedPortfolioValuesTenDays.at(i) - portfolioValue;
//		std::cout << std::setw(20) << std::fixed << (portfolioValue - simulatedPortfolioValuesOneDay.at(i))
//				<< std::setw(20) << (portfolioValue - simulatedPortfolioValuesTenDays.at(i))
//				<< std::endl;
		of << i << ";" << std::fixed << portfolioValueChangesOneDay.at(i)
				<< ";" << portfolioValueChangesTenDays.at(i)
				<< std::endl;
	}
	of.close();
	unsigned VaRIdx = NUMBER_OF_SIMULATIONS / 100;
	std::cout << "MonteCarlo VaR one day: " << std::setprecision(2) << std::fixed << portfolioValueChangesOneDay.at(VaRIdx)
			<< " EUR" << std::endl;
	std::cout << "MonteCarlo VaR ten days: " << std::setprecision(2) << std::fixed << portfolioValueChangesTenDays.at(VaRIdx)
			<< " EUR" << std::endl;

	// calculate arithmetic average of simulated values
	std::vector<double> simAverOneDay(nenEndOfCompanies);
	std::vector<double> simAverTenDays(nenEndOfCompanies);

	for (unsigned n=0; n < NUMBER_OF_SIMULATIONS ; n++)
	{
		for (unsigned c = 0; c < nenEndOfCompanies; c++)
		{
			simAverOneDay.at(c) += mOneDaySimulatedCourseChanges.at(n).at(c);
			simAverTenDays.at(c) += mTenDaySimulatedCourseChanges.at(n).at(c);
		}
	}
	for (unsigned c = 0; c < nenEndOfCompanies; c++)
	{
		simAverOneDay.at(c) /= NUMBER_OF_SIMULATIONS;
		simAverTenDays.at(c) /= NUMBER_OF_SIMULATIONS;
	}
	int width = 20;
	int prec = 6;
	std::cout << std::setw(width) << "Comp." << std::setw(width) << "sim. one day aver."
			<< std::setw(width) << "sim. ten days aver." << std::endl;
	for (unsigned c = 0; c < nenEndOfCompanies; c++)
	{
		std::cout << std::setw(width) << c << std::setw(width) << std::setprecision(prec) << simAverOneDay.at(c)
					<< std::setw(width) << simAverTenDays.at(c) << std::endl;
	}

	// covariance matrices of simulation
	CovarianceMatrix simOneDayCovMatrix;
	CovarianceMatrix simTenDaysCovMatrix;
	double sum,diffI,diffJ;
	for (unsigned i = static_cast<unsigned>(nenDaimler) ; i < static_cast<unsigned>(nenEndOfCompanies); i++) {
		simOneDayCovMatrix.push_back(std::vector<double>(nenEndOfCompanies));
		simTenDaysCovMatrix.push_back(std::vector<double>(nenEndOfCompanies));
	}

	//one day
	for (unsigned i = static_cast<unsigned>(nenDaimler) ; i < static_cast<unsigned>(nenEndOfCompanies); i++) {
		for (unsigned j = static_cast<unsigned>(nenDaimler) ; j < static_cast<unsigned>(nenEndOfCompanies); j++) {
			if (i <= j) {
				sum = 0.0;
				for (unsigned n=0 ; n < NUMBER_OF_SIMULATIONS ; n++) {
					diffI = mOneDaySimulatedCourseChanges.at(n).at(i) - simAverOneDay.at(i);
					diffJ = mOneDaySimulatedCourseChanges.at(n).at(j) - simAverOneDay.at(j);
					sum += (diffI * diffJ) / (NUMBER_OF_SIMULATIONS - 1);
				}
				simOneDayCovMatrix.at(i).at(j) = sum;
			} else {
				simOneDayCovMatrix.at(i).at(j) = simOneDayCovMatrix.at(j).at(i);
			}
		}
	}

	// ten days
	for (unsigned i = static_cast<unsigned>(nenDaimler) ; i < static_cast<unsigned>(nenEndOfCompanies); i++) {
		for (unsigned j = static_cast<unsigned>(nenDaimler) ; j < static_cast<unsigned>(nenEndOfCompanies); j++) {
			if (i <= j) {
				sum = 0.0;
				for (unsigned n=0 ; n < NUMBER_OF_SIMULATIONS ; n++) {
					diffI = mTenDaySimulatedCourseChanges.at(n).at(i) - simAverTenDays.at(i);
					diffJ = mTenDaySimulatedCourseChanges.at(n).at(j) - simAverTenDays.at(j);
					sum += (diffI * diffJ) / (NUMBER_OF_SIMULATIONS - 1);
				}
				simTenDaysCovMatrix.at(i).at(j) = sum;
			} else {
				simTenDaysCovMatrix.at(i).at(j) = simTenDaysCovMatrix.at(j).at(i);
			}
		}
	}

	// dump covariance matrices
	width = 12;
	std::cout << "Covariance matrix for simulation for one days" << std::endl;
	for (unsigned i = static_cast<unsigned>(nenDaimler) ; i < static_cast<unsigned>(nenEndOfCompanies); i++) {
		for (unsigned j = static_cast<unsigned>(nenDaimler) ; j < static_cast<unsigned>(nenEndOfCompanies); j++) {
			std::cout << std::setw(width) << std::setprecision(prec) << simOneDayCovMatrix.at(i).at(j);
		}
		std::cout << std::endl;
	}
	std::cout << "Covariance matrix for simulation for ten days" << std::endl;
	for (unsigned i = static_cast<unsigned>(nenDaimler) ; i < static_cast<unsigned>(nenEndOfCompanies); i++) {
		for (unsigned j = static_cast<unsigned>(nenDaimler) ; j < static_cast<unsigned>(nenEndOfCompanies); j++) {
			std::cout << std::setw(width) << std::setprecision(prec) << simTenDaysCovMatrix.at(i).at(j);
		}
		std::cout << std::endl;
	}
}