Monte Carlo is most useful when you lack analytic tractability or when you have a highly multidimensional problem.
For example, even using simple lognormal and poisson models, there exist path-dependent payoffs or multi-asset computations such that no analytic solution exists and such that any PDE finite difference solution would require 3 or more dimensions. Other times, you are employing a model where the SDE is not solvable, so an apparently one-dimensional problem still ends up forcing you to generate many incremental paths using Euler or Milstein integration.
Cases Where Monte Carlo Is A Poor Idea
- Weakly path-dependent options (e.g. lookbacks): Use PDE or series solutions
- Single-dimensional cases: If your problem is just one dimensional, such as pricing a payoff along the terminal distribution, you should never use Monte Carlo, since numerical quadrature is far superior in this case, even if you just use Riemann sums.
Cases Where Monte Carlo Is A Good Idea
- Strongly path-dependent options such as ratchet range options
- Portfolio risks and exotic baskets where high dimensionality comes into play. CDO tranche protection is a classic example. So are tail risk computations, especially for multi-asset portfolios.
- Intractable models where the terminal distribution cannot be computed, such as some stochastic vol models
To the point about single-dimensional cases -- it sounds like this describes your usage, perhaps because you are using some kind of implied distributional fit to agree with volatility skew. In this case Monte Carlo seems easy, but using a trapezoid rule integrator (or similar) will be as easy and far higher quality by about any measure.
Now Monte Carlo does make it tricky to accurately compute greeks. As with any model, we can compute greeks by using a finite difference "parameter bump", computing our greek
$$
g_\mu =\frac{ V(\dots,\mu+\Delta \mu,\dots) - V(\dots, \mu,\dots)}{\Delta \mu}
$$
but if there is a lot of random noise in those two separate computations of $V()$ then our $g_\mu$ will be inaccurate. Instead it is important to bring the differencing inside the Monte Carlo formula. That is, we don't want to be doing
$$
\hat{g}_\mu =\frac{ \frac1M \sum_{i=1}^M V(x_i,\dots,\mu+\Delta \mu,\dots) - \frac1M \sum_{i=1}^MV(y_i, \dots, \mu,\dots)}{\Delta \mu}
$$
for two separate sample sets $x_i$ and $y_i$. Instead, we want to use the same $x_i$ for both sums, meaning we effectively compute
$$
g_\mu =\frac1{M {\Delta \mu}} \sum_{i=1}^M V(x_i,\dots,\mu+\Delta \mu,\dots) -V(x_i, \dots, \mu,\dots)
$$
and end up with a far more accurate estimate, typically better than our estimate of option value.
I'll make one final note, which is that you feel you "get sound results without relying on the assumptions of the analytical methods". If your terminal distributions are empirically generated, then you are likely to misprice any options because you are not using anything close to a risk-neutral measure. For example, you almost certainly find yourself pricing a forward contract $F$ far higher than the true, arbitrageable value range $F \in [S_0 e^{r_L T}, S_0 e^{r_b T}]$ where $r_b, r_L$ are standard borrow-lend rates.
(Vytautas beat me to some of these points)