Monte Carlo Pricing with Two-State Markov Volatility

# Jump Diffusion Process question


# Jump Diffusion Process question












I have a European call option with time maturity $T=3$ years,$K=50$, and given that $S(t)$ refers to the derivative is being described by the geometric Brownian motion with $S_{0}=100$ and $r = 0.04$.The volatility of the title is stochastic and deviates between two stages $\sigma_-=0.25$ and $\sigma_+ =0.75.$Those two stages can be described as a process of jumps where the time between $0.25 \to 0.75$ is exponentially distributed with $\lambda_+=1$ and the jump $0.75 \to 0.25$ also exponentially distributed with $\lambda_-=3$.So under the risk neutral measure the process is described by $$dS(t) = (r-\sigma_{t}^2/2)S(t)dt+\sigma_{t}S(t)dW(t)$$

My question is how I can model this "two stages" with two different sigmas and different lambdas of exponentially distributed time arrival ?

## Answer by Kermittfrog (score 1, accepted)

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

NB: I reworked the answer to make use of the sojourn time distribution of CTMCs, see below.

#### A sketchy first ansatz

Preliminary remark: Under the proposed setup, the volatility follows a (continuous time) Markov Chain on two volatility states with a high-volatility state ($h$) and a low-volatility state ($l$) and jump intensities $\lambda_{l\to h}$, $\lambda_{h\to l}$. Let us simplify a bit and assume a discrete jump process with time step $\Delta t = \tau / n$ and transition matrix

$$ Q\equiv\begin{pmatrix} P(l_t\to l_{t+\Delta t}) & P(h_t\to l_{t+\Delta t})\\P(l_t\to h_{t+\Delta t}) & P(h_t\to h_{t+\Delta t})\end{pmatrix}\approx \begin{pmatrix} 1-\Delta t\lambda_{l\to h} &\Delta t\lambda_{h\to l}\\\Delta t\lambda_{l\to h} & 1-\Delta t\lambda_{h\to l}\end{pmatrix} $$ Let' simplify the matrix as $$ \begin{pmatrix}P_{ll}&P_{hl}\\P_{lh}&P_{hh}\end{pmatrix} $$

Given some sensible discretization parameter $n$ and some initial state $v_1=h$ or $v_1=l$, we have all ingredients for calculating the distribution of the 'number of times' $\#h, \#l=n-\#h$ the chain is in the high and low volatility states. For example, with $n=5$ and $v_1=l$

$$ \begin{align} P_n(\#l=5)=&P_{ll}^5\\ P_n(\#l=4)=&P_{ll}^4P_{lh}\\ +&4P_{ll}^3P_{lh}P_{hl}\\ P_n(\#l=3)=&P_{ll}^3P_{lh}P_{hh}\\ +&3P_{ll}^2P_{lh}^2P_{hl}\\ +&2P_{ll}^2P_{lh}P_{hl}P_{hh}\\ +&3P_{ll}P_{lh}^2P_{hl}^2\\ +&P_{ll}^2P_{lh}P_{hl}P_{hh}\\ \end{align} $$ ...etc. In theory, we could pre-compute these probabilities for some sensible level $n$ and come up with a (discrete) distribution of $\sigma$ from which we can produce a Merton-like call option formula. But I have the feeling that this will be quite cumbersome, to be honest.

#### A second ansatz

We know that sojourn times of continuous time markov processes are exponentially distributed. For example, given initial state $l$, the first jump time $t_1$ is exponentially distributed with

$$ t_1\sim exp(\lambda_{l\to h}) $$ Then, as above, the probability that we will see at least one jump is $P(t_1<\tau)$, and we will have to calculate some very involved convolutions to come up with the first couple of sensible distributions for the times in each state.

Hence, I suggest a

#### pure MC ansatz .

Fix the usual parameters, especially the time to maturity `tau`, the volatility levels `vol_lo` and `vol_hi`, and some initial state indicator `state0=1` (low) or `state0=2` (high). Then we make use of the exponential distribution of the sojourn times (adjust to your application)

```
set.seed(42)

nSim   <- 1E5
tau    <- 3
K      <- 105
S0     <- 100
rf     <- 0.05
vol_lo <- 0.25
vol_hi <- 0.75
lambda <- c(3,2) # away-from-lo, away-from-hi

sim_time_in_lo <- function(state0){
  t <- 0
  s <- state0
  time_lo <- 0
  while (t<tau){
    dt <- rexp(n=1,lambda[s])
    if ((t+dt)>tau){ dt <- tau - t}
    if (s==1){time_lo <- time_lo + dt }
    if (s==1){s<-2} else {s <-1}
    t <- t+ dt
  }
  time_lo
}
tau_lo     <- sapply(1:nSim,function(i){sim_time_in_lo(1)})
tau_hi     <- tau - tau_lo
total_var  <- tau_lo * vol_lo^2 + tau_hi * vol_hi^2
drift      <- rf * tau-0.5*total_var
randomness <- sqrt(total_var)*rnorm(nSim,)

mean(exp(-rf*tau)*pmax(S0*exp(drift + randomness)-K,0))
```

with output:

```
[1] 41.58493
```

I hope that helped a bit? I'll put in some more thought during the next days; and maybe somebody else has an idea.