Let's assume I have an arbitrary option that I can price using Monte-Carlo simulation. What is the general approach (i.e. without relying on specific option type) to calculating the greeks in this case?
Edit: I woud like to add a few links on the topic that I found useful:
You need to compute your greeks as finite differences, but the full procedure may be pretty tricky. I will use vega $\aleph$ as the example here. Let's begin by designating your Monte Carlo estimator as a function $V(\sigma,s,M)$ where $\sigma$ is the volatility as usual, $s$ is the seed to your random number generator, and $M$ is the sample count.
To begin with, recall that the Monte Carlo estimate of any value converges with the square root of the sample count. In particular, if you choose, say, $M=100$, you can run your estimator $N=500$ times to get estimate $\{V_n\}_{i=1}^{500}$, obtaining the standard deviation $\Sigma_{100}$ of those estimates.
Having done this, we now know the standard error of the estimator for any $M$ to be
$$ e_M \approx \Sigma_{100} \sqrt{\frac{100}{M}} $$
There are three possible cases:
In the first case, you can use the fact that $s$ has been controlled to get a reasonable estimate of vega with relatively little extra work.
Find an $M$ such that the error in option price $e_M$ is tolerable. Choose a seed $s_0$ and a small increment $\Delta\sigma$ in the volatility, and compute
$$ \aleph^{(1)} = \frac{V(\sigma+\Delta\sigma,s_0,M)-V(\sigma-\Delta\sigma,s_0,M)}{2\Delta\sigma} $$ and $$ \aleph^{(2)} = \frac{V(\sigma+\frac12\Delta\sigma,s_0,M) - V(\sigma-\frac12\Delta\sigma,s_0,M)}{\Delta\sigma} $$
If $\aleph^{(1)} \approx \aleph^{(2)}$ then you have a good estimate and you are done.
The reason this works so nicely is that, by controlling the seed, our difference computations
$$ \delta=V(\sigma+\Delta\sigma,s_0,M)-V(\sigma-\Delta\sigma,s_0,M) $$ are direct Monte Carlo estimators of the vega, since the shared seed implies the samples $x_i$ match in the difference of sums. That is $$ \delta = ( \frac1M \sum_{i=1}^M f(x_i, \sigma+\Delta\sigma) ) -( \frac1M \sum_{i=1}^M f(x_i, \sigma-\Delta\sigma) ) \\ =\frac1M \sum_{i=1}^M f(x_i, \sigma+\Delta\sigma)-f(x_i, \sigma-\Delta\sigma) $$
The second case where you cannot control the seed, on the other hand, is rather more difficult. Here, you will have a different error $e$ to the true value every time you run the function.
For brevity, let's let $$ e_\pm = V(\sigma\pm \Delta\sigma,s_\pm,M). $$
Of course we do not know the value of $e_\pm$ or of $s_\pm$, but we do at least have our estimate of the size of $e_\pm$ as noted above. Therefore, the error in $\delta$ is approximately $e_M \sqrt{2}$. You need to choose $M$ so large that $$ \delta \gg e_M \sqrt{2}. $$
Not knowing the value of $\delta$ a priori makes this difficult, but usually in a trading context one can specify an acceptable absolute error $\epsilon$ in vega. In that case, we can demand $$ \epsilon < \frac{e_M \sqrt{2}}{\Delta\sigma} $$ which translates to $$ M > \Sigma_{100}^2 {\frac{200}{\epsilon^2 \Delta\sigma^2}}. $$
The third case, where you can control neither the random seed $s$ nor the sample count $M$ should be treated as the second case above. You simply treat each run of the algorithm as a single sample.
The most general answer is to shift your input to approximate the first derivative. Given that you need Monte Carlo to price this, it may get expensive. But that's the way it goes as when you have no analytical solutions as there aint't no free lunch ...
Download the code from http://fmsoption.codeplex.com to see how to do that for vanilla options. You are right, you need implementations for transcendental functions that are written for dual numbers. You will find them in the fmsdual project.
If you just want browse some source code, see http://fmsoption.codeplex.com/SourceControl/changeset/view/10924#145366. Note that eps is machine epsilon ~= 2e-16. (!)
if the pay-off is continuous, the standard approach is to use the path-wise method also known as IPA. This essentially means that you differentiate along each path. It is the limit as the bump size goes to zero of finite differencing.
The main downside of this method is that the differentiation can be fiddly and slow. The Smoking adjoints paper you mention made the observation that using adjoint/automatic differentiation makes it fast and indeed there are packages that will do the differentiating automatically for you. There is a survey article by Homescu on this topic.
Adjoints and Automatic (Algorithmic) Differentiation in Computational Finance
Cristian Homescu
I have also written too many papers on this. See my SSRN page http://papers.ssrn.com/sol3/cf_dev/AbsByAuth.cfm?per_id=550354
I also devote a chapter of More Mathematical Finance to this.
When the pay-off is discontinuous or you want Hessians life is more complicated. Likelihood ratio is one approach but can lead to large variances. This is particularly the case for vegas.
See Glasserman (2003) Monte Carlo Methods in Financial Engineering for general discussion.
I have done various papers on the discontinuous case and on Hessians. My current favourites are
http://ssrn.com/abstract=2431580
http://ssrn.com/abstract=2011690
Essentially, these work by doing change of variables that remove the discontinuities and so then the pathwise method can be applied. They are designed to yield the most pathwise mix of likelihood ratio and pathwise that avoids the discontinuity.
Another way to do this is to use dual numbers. http://fmsdual.codeplex.com.
They let you calculate an arbitrary number of derivatives while running a single Monte Carlo. Here is an example of how to use it:
// Monte Carlo derivatives
void fms_test_monte(size_t N)
{
::srand(static_cast<unsigned int>(::time(0)));
double a = 0.5;
dual::number<double,3> A(a, 1);
dual::number<double,3> E(0.,1);
for (int i = 0; i < N; ++i) {
double x = 1.0*rand()/RAND_MAX;
E = E + (x - A)*(x - A);
}
E = E/(1.*N);
// X uniform [0,1]
// E(X - a)^2 = 1/3 - 2a 1/2 + a^2
ensure (fabs(E._(0) - (1./3 - a + a*a)) < sqrt(1./N));
// d/da E(X - a)^2 = -2 E(X - a) = 2a - 1
ensure (fabs(E._(1) - (2*a - 1)) < sqrt(1./N));
// d^2/da^2 E(X - a)^2 = 2
ensure (fabs(E._(2) - 2) < sqrt(1./N));
}