EDIT 2: I found the problem(s) and the prices seem to behave as expected now. For anyone interested there was a bug when normalizing the dependant ranom normal variates used in the simulation, so while they had the correct correlation one of them had a standard deviation of 1 and the other a standard deviation much greater than 1. Causing the price to not drop (and even increase) even as correlation increased. The parity relation suggested by @ir7 seems to hold now which gives me confident that all is well.
I'm valuing a rainbow option numerically with a monte-carlo simulation and I'm getting some unexpected results. The price of a best-of call option is decreasing in price as correlation up to a certain point where it starts to increase, contrary to my intuition. The worst-of option case is much more well behaved with being an increasing function of the correlation as expected. Since I'm valuing them in much the same way (just taking min(...) instead of max(...) in my code) I'm very confused as to what could be wrong, or is it possible for the price to behave in this way? If it's completely unreasonable, does anyone want to hazard a guess as to why my calculations might be breaking down as the correlation increases? The program is written in C++ so if anyone fluent in C++ wants to have a look at my code for something erroneous I'd be more than happy to post it.
EDIT 1: After troubleshooting a bit with the help of @ir7 it seems there's something off with my Monte Carlo simulation for several assets (the single asset case works fine). Posting (some of) my C++ code below for the simpler problem of valuing an out-performance option which as a closed form solution for him (and anyone else who wants to of course) to check out and give help. The caclulations again seem to break down once the correlation goes above ~0.5, see below image. If there are some function calls used that you want me to explain or post the code for I'd be happy to do so, for now I'll try to keep it somewhat bare:
The class and function that does the actual valuation:
MonteCarloOutPerformanceOptionFunction::MonteCarloOutPerformanceOptionFunction(std::string uniqueIdentifier_, int nominal_, std::vector<double> S0_vect, std::vector<Wrapper<PayOff>> ThePayOffVect_, double r_, std::vector<double> d_vect_, std::vector<double> impvol_vect_, std::vector<std::vector<double>> covMatrix_, double TTM_, unsigned long numberOfPaths_)
: r(r_), S_vect(S0_vect), ThePayOffVect(ThePayOffVect_), d_vect(d_vect_), covMatrix(covMatrix_), valuationFunction(uniqueIdentifier_, TTM_, nominal_), numberOfPaths(numberOfPaths_), impvol_vect(impvol_vect_)
{
if (covMatrix.size() != S_vect.size())
throw("Missmatched Covariance matrix and initial spot values array sizes in OutPerformance Option");
if (2 != S_vect.size())
throw("More than two equities specified in OutPerformance Option");
}
void MonteCarloOutPerformanceOptionFunction::ValueInstrument()
{
std::vector<MJArray> correlatedNormVariates = GetArraysOfCorrelatedGauassiansByBoxMuller(numberOfPaths, covMatrix);
std::vector<StatisticAllPaths> thesePathGatherers;
for (unsigned long i = 0; i < S_vect.size(); i++)
{
StandardExcerciseOption thisOption(ThePayOffVect[i], TTM);
StatisticAllPaths onePathGatherer;
thesePathGatherers.push_back(onePathGatherer);
OneStepMonteCarloValuation(thisOption, S_vect[i], impvol_vect[i], r, d_vect[i], numberOfPaths, correlatedNormVariates[i], thesePathGatherers[i]);
}
f = 0;
for (unsigned long i = 0; i < numberOfPaths; i++)
{
std::vector<double> outcomes;
outcomes.reserve(S_vect.size());
for (unsigned long j = 0; j < S_vect.size(); j++)
{
outcomes.push_back(thesePathGatherers[j].GetOneValueFromResultsSoFar(i));
}
f += std::max(outcomes[0] - outcomes[1], 0.0);
}
f *= ((double)nominal / numberOfPaths);
return;
}
The Monte Carlo simulation function being called at OneStepMonteCarloValuation (this seems to work fine for single asset options like vanilla calls/puts)
void OneStepMonteCarloValuation(const StandardExcerciseOption& TheOption, double Spot, double Vol, double r, double d, unsigned long NumberOfPaths, MJArray normVariates, StatisticsMC& gatherer)
{
if (normVariates.size() != NumberOfPaths)
throw("mismatched number of paths and normal variates");
//Pre-calculate as much as possible
double Expiry = TheOption.GetExpiry();
double variance = Vol * Vol * Expiry;
double rootVariance = sqrt(variance);
double itoCorrection = -0.5 * variance;
double movedSpot = Spot * exp((r-d) * Expiry + itoCorrection);
double thisSpot;
double discounting = exp(-r * Expiry);
for (unsigned long i = 0; i < NumberOfPaths; i++)
{
thisSpot = movedSpot * exp(rootVariance * normVariates[i]);
double thisPayoff = TheOption.OptionPayOff(thisSpot);
gatherer.DumpOneResult(discounting * thisPayoff);
}
return;
}
The StatisticAllPaths class which is used as input in the simulation that collects all the final values of the simulation
StatisticAllPaths::StatisticAllPaths(const unsigned long minimumNumberOfPaths) : PathsDone(0)
{
ResultList.reserve(minimumNumberOfPaths);
}
void StatisticAllPaths::DumpOneResult(double result)
{
ResultList.push_back(result);
PathsDone++;
}
const double& StatisticAllPaths::GetOneValueFromResultsSoFar(unsigned long index) const
{
return ResultList[index];
}
The PayOffVect used is used here to take the payoff of each path in the MC valuation function, but since we're just collecting all the paths here and processing them later (in the last part of the main valuation class) it doesn't really do anything here. It it used in this case just to make the outperformance relative values with this inherited class:
PayOffRelPerformance::PayOffRelPerformance(double startValue_) : startValue(startValue_)
{
}
double PayOffRelPerformance::operator()(double spot) const
{
return spot / startValue;
}
The GetArraysOfCorrelatedGauassiansByBoxMuller does the job of generating the vectors of normal variates that will be used in the simulation. I have checked that the Cholezky Matrix is correct for real cases, and I have also checked that the outputted normal variates are in fact dependent with the correlation implied by the covariance Matrix.
std::vector<MJArray> GetArraysOfCorrelatedGauassiansByBoxMuller(unsigned long numberOfVariates, std::vector<std::vector<double>> covMatrix)
{
//Calculate the cholezky Matrix
std::vector<std::vector<double>> cholezkyMatrix = Cholesky_Decomposition(covMatrix);
//Fix the size of the arrays to contain correlated normal variates
std::vector<MJArray> corrNormVariatesVector(cholezkyMatrix.size());
for (unsigned long j = 0; j < corrNormVariatesVector.size(); j++) {
corrNormVariatesVector[j].resize(numberOfVariates);
corrNormVariatesVector[j] = 0;
}
//calculate correlated normal variates and fill the arrays with values
MJArray NormVariates(cholezkyMatrix.size());
for (unsigned long k = 0; k < numberOfVariates; k++) {
for (unsigned long i = 0; i < cholezkyMatrix.size(); i++)
{
NormVariates[i] = GetOneGaussianByBoxMuller();
for (unsigned long j = 0; j < cholezkyMatrix[i].size(); j++) {
corrNormVariatesVector[i][k] += cholezkyMatrix[i][j] * NormVariates[j];
}
corrNormVariatesVector[i][k] /= cholezkyMatrix[i][i]; //normalize the random variates
}
}
return corrNormVariatesVector;
}