# Simulating the joint dynamics of a stock and an option

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?

-
Just to be clear, are you assuming there exists an options pricing function BSM(S,t,{params}) or are options prices also supposed to come from the Monte Carlo simulation? –  Brian B Feb 29 '12 at 20:12
I've updated the question, I assume the existence of the options pricing function $BSM(S,t,{params})$. –  Bob Jansen Feb 29 '12 at 20:50

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.$$
Thanks, I noticed the $T - \tau$ and wondered whether it was correct under a different interpretation of the variables. Nice to know it's not ;) I've implemented most of it over the weekend and will show the result (appendix in my thesis) and MATLAB code soon, I hope. You've really been a great help, quadrature is new terrain for me and your answer was a great start. –  Bob Jansen Mar 5 '12 at 15:46