How to get greeks using Monte-Carlo for arbitrary option?

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:

Answer

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:

1. You can control the random seed $s$, or the set of random samples, used by the Monte Carlo estimator


2. You cannot control $s$.


3. You cannot even control the sample count $M$.

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.