Simulating Dividend-Paying Stocks with the Milstein Scheme

# MonteCarlo simulation of stock prices using milstein scheme with dividend yield?


# MonteCarlo simulation of stock prices using milstein scheme with dividend yield?












While performing a montecarlo simulation of stock prices using the milstein scheme is it possible to take into account the dividend yield into the simulation itself somehow, if we are given a continuous dividend yield?

Or is this something that has to be considered while valuing various derivatives using those simulated paths?

Am fairly new to quant finance, so pardon me if this turns out to be a stupid question.

## Answer by user16891 (score 3)

https://quant.stackexchange.com/a/19003

It is straightforward to include dividends into the model if it can be assumed that the dividend payment is a continuous dividend yield, $q$. Under $Q$ measure, in the Black-Scholes Model, Heston Model etc, $r$ is replaced by $r − q$., Here, We are going to simulate underlying asset in the Black-Scholes model by Milstein Method. Indeed, we assume underlying asset follows the Ito process as described by the following stochastic differential equation $$dS_t=(r-q)S_t dt+\sigma\, S_tdW_t$$ Apply Milstein discretization to this equation,we have $$S_{t+\Delta t} = S_t+(r-q)S_t \Delta t+\sigma\,S_t \sqrt{\Delta t}\,Z+\frac{1}{2}\sigma^2\Delta t(Z^2-1)$$ where $Z$ is distributed as standard normal. We then retain the last stock price from each stock price path and obtain the payoff of the European option at expiry, take the average over all stock price paths and discount back to time zero. Hence, for example, the call $C(K)$ and $C(K,S_T,r,\sigma)$ is $$C(K,S_T,r,\sigma)=e^{-rt}\frac{1}{N}\sum_{i=1}^{N}max\{S_T^{(i)}-K, 0\}$$

## Answer by Petr Andreev (score 3)

https://quant.stackexchange.com/a/64227

The right formula: $$ S_{t+\Delta t} = S_t+(r-q)S_t \Delta t+\sigma\,S_t \sqrt{\Delta t}\,Z+\frac{1}{2}\sigma^2\Delta t(Z^2-1)*S_t $$

We can extract the formula from the Brownian motion equation (Wiener process): $$ dS_t = (r-q) S_t dt + \sigma S_t dW_t $$ where $ W_t $ is Wienere process.

Applying Itô's lemma (r, q and $\sigma$ are constants) with $ f(S_t) = \ln(S_t) $ gives:

$$ df(S_t) = f^\prime(S_t)*dS_t + \frac{1}{2}f^{\prime\prime} (S_t) *(dS_t)^2 \\ $$ where $ (dS_t)^2 $ is the quadratic variation of the SDE: $ (dS_t)^2 = \, \sigma^2 \, S_t^2 \, d W_t^2 + 2 \sigma S_t^2 (r-q) \, d W_t \, d t + (r-q)^2 S_t^2 \, d t^2 $

When $ dt \to 0 $, $ dt $ converges to 0 faster than $ dW_t $, since $ dW_t^2 = O(dt) $.

So the above infinitesimal can be simplified by: $ (dS_t)^2 = \sigma^2 S_t^2 dt $

$ df(S_t) = \frac{1}{S_t}\,dS_t + \frac{1}{2} (-S_t^{-2}) (S_t^2\sigma^2\,dt) = \frac{1}{S_t} \left( (r-q) S_t\,dt + \sigma S_t\,dW_t \right) - \frac{1}{2}\sigma^2\,dt = \\ =\left (r-q-\frac{1}{2}\sigma^2\ \right )\,dt + \sigma dW_t $

$$ d(ln S_t) = (r-q - \frac{1}{2}\sigma^2) dt + \sigma dW_t $$ $$ S_{t+\Delta t} = S_t * e^{\int\limits_t^{t+\Delta t}{(r - q - \frac{1}{2}\sigma^2)}dt + \int\limits_t^{t+\Delta t}{\sigma dW_u}} $$ taking integrals and decomposing the exponent to Maclaurin series up to the third term $$ ≈ S_t * [1 + (r-q)\Delta t + \sigma \Delta W_t + \frac{1}{2}\sigma^2((\Delta W_t)^2 - \Delta t)] $$ where $ \Delta W_t = Z(0,1)*\sqrt{\Delta t}$, where $ Z(0,1) $ is a random variable from normal distribution $$ S_{t+\Delta t} = S_t+(r-q)S_t \Delta t+\sigma\,S_t Z \sqrt{\Delta t} + \frac{1}{2}\sigma^2(Z^2-1) \Delta t*S_t $$

p.s.: user16891 forgot to multiply the last term by $ S_t $

p.s.s.: using the equation $ S_{t+\Delta t} = S_t * e^{(r - q) \Delta t + \sigma Z \sqrt{\Delta t} + \frac{1}{2}\sigma^2 (Z^2-1) \Delta t} $ for Monte-Carlo simulation can be faster (in python).

p.s.s.s.: more theory about $ dW $ you can find here (Course Notes on SDEs from Moscow State University, Russian language)