Simulating the joint dynamics of a stock and an option

Question

I want to know the joint dynamics of a stock and it's option for a finite number of moments between now and $T$ the expiration date of the option for a number of possible paths.

Let $r_{\mathrm{s}}$ and $r_{\mathrm{o}}$ denote the return on the stock and the option. Then I'm interested in knowing $\mathrm{E}_t([r_{\mathrm{s}}; r_{\mathrm{o}}])$ the expectation of the return from $t$ to $t+1$ and $\mathrm{Var}_t(r_{\mathrm{s}}, r_{\mathrm{o}})$ the variance from $t$ to $t+1$. The expected return and volatility of the stock are known.

My first idea is to use Monte Carlo with the following pseudocode:

N <- number of paths
T <- number of moments
M <- number of subpaths
S <- current stock price
for i = 1 to N:
    S_0 <- S
    for t = 0 to T-1:
        for j = 1 to M:
            S_{t+1,j} = f(S_t)
            O_{t+1,j} = BSM(S_{t,j})
        S_{t+1} <- mean(S_{t+1,j})
        E_{i,j} <- mean(S_{t+1,j}, O_{t+1,j})
        V_{i,j} <- var(S_{t+1,j}, O_{t+1,j})
return (E, V)

where $S_t$ is the current stock price, $f(S_t)$ gives a realization of stock price at the next moment given the current stock price, let's assume geometric Brownian Motion, $\textrm{BSM}(S_t)$ gives the option price given the current stock price and some arbitrary parameters using the BSM formula and E and V the values I'm interested in.

This can probably be optimized by a number of clever ideas such as reusing the draws from the probability distribution and discretizing the state space and use a memoized BSM. This is, however, not what I'm looking for. I rather calculate the mean and variance directly, the question is: how?

Answer 1

In this scenario, the "joint dynamics" are trivially computed since the option value is a known deterministic function of stock price. For example, the mean of the option value for time $\tau$ is $$ \mu_O = \int_0^\infty BSM( S_\tau ) p(S_\tau) dS_\tau $$ which is best computed using quadrature as available in standard numerical libraries like scipy. The function $p(S_\tau)$ would typically be the Black-Scholes probability density $$ \frac{n( d_2(S_0,S_\tau) )} {S_\tau \sigma \sqrt{\tau} }. $$ with $S_\tau$ taking the place of strike $K$ in the formula for $d_2()$.

Similarly, the variance of the option value for time $\tau$ is $$ \int_0^\infty (BSM( S_\tau ) - \mu_O)^2 p(S_\tau) dS_\tau $$ and covariance of option value with stock price would be $$ \int_0^\infty (BSM( S_\tau ) - \mu_O)(S_\tau-\mu_S) p(S_\tau) dS_\tau. $$