# 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:

-
That vibrato Monte Carlo thesis is very interesting. –  Brian B Apr 27 '12 at 15:00

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.

-

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 ...

-
Numerical derivatives are iffy business, but I agree that it seems to be your best choice. As you probably know; be aware of the how the precision decreases quickly(!) as higher orders are measures. –  AdAbsurdum Apr 24 '12 at 7:34

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. (!)

-
For the curious, the "dual number" approach referenced here is an automatic differentiation package, where the bookkeeping mostly handled by C++ templates and extensions to the standard numeric types. Derivatives of transcendental functions are handled by automatic differentiation of the numerical analytic series approximations used to calculate function values. (Please correct any mistakes I have made in that) –  Brian B May 10 '12 at 14:12
I beleive AD refers to techniques for automatically generating functions for the derivatives. Dual numbers don't do that for you. –  Keith A. Lewis May 11 '12 at 14:05

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));
}

-
A novel method. It appears mathematically equivalent to finite differences. Also worth noting is that all samples need to remain in memory, at least if I read this stuff right. It would be informative to see code applied to a nontrivial estimator function. –  Brian B Apr 27 '12 at 15:36
You are reading it wrong. Not all samples need to remain in memory. It is not mathematically equivalent to finite differences, you seem to be completely missing the point of dual numbers. You have the source code, let me know if you need help applying it to a non-trivial estimator. Monte Carlo aside, dual numbers allow you to calculate derivatives down to machine precision. –  Keith A. Lewis May 7 '12 at 9:56
You're right, I now see it is not finite differences. I would find the simplicity of your example far more convincing if you demonstrated calculating, say, delta of an average strike option under the the Black-Scholes stochastic model. It seems to me the transcendental functions involved make this difficult even for vanilla options. Path dependencies will make the problem much worse. –  Brian B May 7 '12 at 16:24
I'll also note that, as @Dirk and I read the question, Alexey does not necessarily have the source code to the option pricer. –  Brian B May 7 '12 at 16:25