Monte Carlo Greeks for Fixed Strike Asian Call

I am interested in pricing an European-style fixed strike asian call with payoff $$\max(A(S)-K;0)$$, where $$A(S)=\frac{1}{n}\sum_{i=1}^nS(t_i)$$ is a discrete arithmetic average and $$K$$ is the strike price.

Assuming an arbitrage-free and complete market, the fundamental theorem of asset pricing tells us that the arbitrage-free price at time $$t=0$$ is given by: $$V(0)=E^{\mathbb{Q}}(\max(A(S)-K) \vert {\cal F}_0 )$$ I have no idea whether there exists an analytic solution or not, so I decided to use MC by implementing the following pseudo-code in python (I omit the code so as not to prolong the question).

Under the Black-Scholes assumptions, let $$m$$ be the number of paths, $$n$$ be the number of intervals per path and $$\delta t= \frac{T}{n}$$, then:

1. Simulate geometric brownian motion under $$\mathbb{Q}$$ measure $$S_i(t+1)=S_i(t) \exp \left(\left(r-\frac{\sigma^2}{2}\right)\delta t+\sigma \sqrt{\delta t}Z_t\right)$$ where $$Z_t \sim{\cal N}(0,1)$$ for $$i \in [1,m]$$ and $$t \in [0,n]$$.
2. Calculate option payoff $$X_i=\max(A_i(S_i)-K)$$ and set $$V_i(0,K,T,\sigma,r,S(0))=e^{-rT}X_i$$
3. Calculate sample average $$V(0,K,T,\sigma,r,S(0))=\frac{1}{m}\sum_{i=1}^m V_i(0,K,T,\sigma,r,S(0))$$

However, when I want to calculate higher order greeks, especially those involving mixed derivatives, things get messy because I have to repeat this process several times and change the parameters slightly. For instance, calculating DdeltaDvol using finite differences yields:

\begin{align*} DdeltaDvol &= \frac{1}{4 \Delta S \Delta \sigma} [V(0,K,T,\sigma+\Delta \sigma,r,S(0)+\Delta S)-V(0,K,T,\sigma-\Delta \sigma,r,S(0)+\Delta S) \\ &\quad -V(0,K,T,\sigma+\Delta \sigma,r,S(0)-\Delta S)+V(0,K,T,\sigma-\Delta \sigma,r,S(0)-\Delta S)] \end{align*} The fact that I have to simulate the whole asset path thousands of times and a with small $$\delta t$$ makes the whole approach computationally intensive.

Does anyone know an alternative approach which is not so computationally intensive ?

• Yes, I forgot to mention that. I fix the seed and then re use the variabels. But I don’t see an easy way to calculate $V(0,\sigma + \Delta \sigma,S(0)+\Delta S)$ from $V(0,\sigma, S(0))$ in a vectorized version, so I have to re-run the loops.