We 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)\]
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.
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 to 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)
nSim <- c(1000,10000,100000)
tau <- 3
K <- c(50,100,150)
S0 <- 100
rf <- 0.04
vol_lo <- 0.25
vol_hi <- 0.75
lambda <- c(3,1) # 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
}
##------------ for k = 50 -------------------
tau_lo <- sapply(1:nSim[1],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[1],)
V <- exp(-rf*tau)*pmax(S0*exp(drift + randomness)-K[1],0)
alpha <- qnorm(1-.05/2)
## Evaluate the European call value
AV <- mean(V) # mean of call price
AVdev <- alpha * sd(V) / sqrt(nSim[1]) # standard error
AV_Upper <- AV + AVdev
AV_lower <- AV - AVdev
Jump_sim1 <- c(AV_Upper,AV,AV_lower)
Jump_sim1## [1] 69.43825 61.47526 53.51226