[R-SIG-Finance] CEV Model & MC simulation

[Francesco Bianchi]

Hi Everyone,

here below the assignment I should solve:


We assume that the underlying process ST follows a CEV model.

  *   Evaluate the Expected Value of the following payoff using the MC method. Write an efficient code. (Hint: system.time() r function could be useful).

max{ST - K, 0} 1 a_1 ≤ ST_1 ≤ b_1, 1 a_1 ≤ ST_2 ≤ b_1, 1 a_1 ≤ ST_3 ≤ b_1

where a_1 < b_1 and 1 is an indicator function, T1 =0.25*T, T2 =0.5*T, T3 =0.75*T

(Hint: The final time T in your code must be an input fixed by the user. If the user chooses T = 10 than T1 =10*0.25 = 2.5 and so on)

  *   Compute the upper and lower bound of the MC price
  *   For each model parameter write a code for evaluating the corresponding sensitivity measure

Below the code used. It sounds too advanced for the course I am enrolled in (it was produced based on advanced reading).
Any way to make it more row and student friendly?

At the end I introduced the sensitivity analysis, what about developing it through a graph analysis instead?

Thanks,

Francesco
_______

#Clean the global environment
rm(list=ls())

#set seed: important to fix when sensitivity analyses are performed
set.seed(123)

cev = function (N,T,x0,gamma,sigma,mu,alpha1,beta1,K,risk_free) {

  #Initialization of variables
  delta.xt=matrix(NA,T,1)     #Vector of dx_t
  xt=matrix(NA,T+1,1)         #Vector of x_t
  xt.condition=matrix(NA,4,N) #Matrix of x_T1, x_T2, x_T3 & x_T (for each scenario - in columns), input of the payoff function
  x.indicator=matrix(NA,3,N)  #Matrix of indicator function (rows of the same column stand for the different conditions in the same scenario)
  payoff=matrix(NA,1,N)       #Vector of payoff, results of the MC simulation

  #Setting initial value of the process, the same for each simulation
  xt[1,1]=x0

  #Definition of T_j condition
  #Ceiling is the rounding function the nearest higher integer
  T1=ceiling(0.25*T)
  T2=ceiling(0.50*T)
  T3=ceiling(0.75*T)

  #Stop the process if alpha1 is higher than beta1
  if (alpha1>beta1) {

    print("Error message: Alpha1 is higher than Beta1 -> not feasible condition -> the MC is stopped without results")

  } else {

    #####Generate the CEV stochastic process
    #Cycle on number of scenario
    for (j in 1:N) {

      #Cycle on time
      for (i in 1:T) {

        #Calculate dx_t and x_t
        #Remember that x_t is 1 unit longer than dx_t as x_t also includes x0 in first position
        #dt is fixed equal to 1
        #J is not used in the following equations, as saving all the dx_t and x_t in each scenario is cumbersome
        delta.xt[i,1]=mu*xt[i,1]*1+sigma*(xt[i,1]^gamma)*rnorm(1,0,1)
        xt[i+1,1]=xt[i,1]+delta.xt[i,1]

      }

      #Save the values of the process in the 4 time units
      xt.condition[1,j]=xt[T1+1,1]
      xt.condition[2,j]=xt[T2+1,1]
      xt.condition[3,j]=xt[T3+1,1]
      xt.condition[4,j]=xt[T+1,1]

      #Print the current number of simulation to monitor the status of MC simulation
      print(paste("Simulation no.",j,"of",N,"completed",sep=" "))

    }

    #Perform the indicator function condition (TRUE=1, FALSE=0)
    x.indicator[1,]=((xt.condition[1,]>=alpha1) || (xt.condition[1,]<=beta1))
    x.indicator[2,]=((xt.condition[2,]>=alpha1) || (xt.condition[2,]<=beta1))
    x.indicator[3,]=((xt.condition[3,]>=alpha1) || (xt.condition[3,]<=beta1))

    #Cycle on number of simulations in order to perform the distribution of payoff
    for (j in 1:N) {

      payoff[1,j]=max(0,xt.condition[4,j]-K)*x.indicator[1,j]*x.indicator[2,j]*x.indicator[3,j]

    }

    #Calculate expected value of the option and its bounds at 95% confidence interval
    expected.value.MC=mean(exp(-risk_free*T/250)*payoff[1,])
    lower.MC = expected.value.MC-1.96*sqrt(var(exp(-risk_free*T/250)*payoff[1,]))/sqrt(N)
    upper.MC = expected.value.MC+1.96*sqrt(var(exp(-risk_free*T/250)*payoff[1,]))/sqrt(N)

  }

  results=list(Expected.Value=expected.value.MC,
               Lower.Bound=lower.MC,
               Upper.Bound=upper.MC)

  return(results)

}

######MonteCarlo speed of convergence

#Reference Simulation

benchmark=cev(N=10000,T=200,x0=2,gamma=1.1,sigma=0.01,mu=0.00001,alpha1=1.4,beta1=2.6,K=2,risk_free=0.01)
benchmark$Expected.Value
benchmark$Lower.Bound
benchmark$Upper.Bound
#In questi tre risultati che stamper? a schermo sono il prezzo di riferimento, il lower bound del MC e l'upper bound del MC.
#L'intervallo di confidenza scelto per il MC ? al 95%

#molto lungo il prossimo codice, si può togliere.
m <- c(10, 50, 100, 150, 250, 350, 500)
p1 <- NULL
err <- NULL
nM <- length(m)
repl <- 100
mat <- matrix(, repl, nM)
for (k in 1:nM) {
  tmp <- numeric(repl)
  for (i in 1:repl) tmp[i] <- cev(N = m[k],T=200,x0=2,gamma=1.1,sigma=0.01,mu=0.00001,alpha1=1.4,beta1=2.5,K=2,risk_free=0.01)$Expected.Value
  mat[, k] <- tmp
  p1 <- c(p1, mean(tmp))
  err <- c(err, sd(tmp))
}
colnames(mat) <- m

minP <- min(p1 - err)
maxP <- max(p1 + err)

plot(m, p1, type = "n", ylim = c(minP, maxP), axes = F, ylab = "MC price", xlab ="MC replications")
lines(m, p1 + err, col = "blue")
lines(m, p1 - err, col = "blue")
axis(2, benchmark$Expected.Value, "Reference price")
axis(1, m)
boxplot(mat, add = TRUE, at = m, boxwex = 15, col = "orange", axes = F)
points(m, p1, col = "blue", lwd = 3, lty = 3)
abline(h = benchmark$Expected.Value, lty = 2, col = "red", lwd = 3)


#Sensitivity Analysis
#questo è un po' lungo ma si puo' togliere

#Greeks with finite difference method

#Delta on X0
delta_price=0.02
looping=100
price_up = matrix(, looping, 1)
price_down = matrix(, looping, 1)

price_stable= matrix(,looping, 1)

for (k in 1:looping) {
  result=cev(N=100,T=200,x0=2,gamma=1.1,sigma=0.01,mu=0.00001,alpha1=1.4,beta1=2.6,K=2,risk_free=0.01)$Expected.Value
  price_stable[k,1]=result
}

for (k in 1:looping) {
  result=cev(N=100,T=200,x0=2+delta_price,gamma=1.1,sigma=0.01,mu=0.00001,alpha1=1.4,beta1=2.6,K=2,risk_free=0.01)$Expected.Value
  price_up[k,1]=result
}

for (k in 1:looping) {
  result=cev(N=100,T=200,x0=2-delta_price,gamma=1.1,sigma=0.01,mu=0.00001,alpha1=1.4,beta1=2.6,K=2,risk_free=0.01)$Expected.Value
  price_down[k,1]=result
}

#Central derivatives
delta_greek_x0=(mean(price_up[,1])-mean(price_down[,1]))/(2*delta_price)

#Vega
delta_vol=0.001
looping=100
price_up = matrix(, looping, 1)
price_down = matrix(, looping, 1)

for (k in 1:looping) {
  result=cev(N=100,T=200,x0=2,gamma=1.1,sigma=0.01+delta_vol,mu=0.00001,alpha1=1.4,beta1=2.6,K=2,risk_free=0.01)$Expected.Value
  price_up[k,1]=result
}

for (k in 1:looping) {
  result=cev(N=100,T=200,x0=2,gamma=1.1,sigma=0.01-delta_vol,mu=0.00001,alpha1=1.4,beta1=2.6,K=2,risk_free=0.01)$Expected.Value
  price_down[k,1]=result
}

vega_greek=(mean(price_up[,1])-mean(price_down[,1]))/(2*delta_vol)

#Gamma_X0 (derivata seconda del sottostante)
delta_price=0.02
looping=100
price_stable= matrix(, looping, 1)
price_up = matrix(, looping, 1)
price_down = matrix(, looping, 1)

for (k in 1:looping) {
  result=cev(N=100,T=200,x0=2,gamma=1.1,sigma=0.01,mu=0.00001,alpha1=1.4,beta1=2.6,K=2,risk_free=0.01)$Expected.Value
  price_stable[k,1]=result
}

for (k in 1:looping) {
  result=cev(N=100,T=200,x0=2+delta_price,gamma=1.1,sigma=0.01,mu=0.00001,alpha1=1.4,beta1=2.6,K=2,risk_free=0.01)$Expected.Value
  price_up[k,1]=result
}

for (k in 1:looping) {
  result=cev(N=100,T=200,x0=2-delta_price,gamma=1.1,sigma=0.01,mu=0.00001,alpha1=1.4,beta1=2.6,K=2,risk_free=0.01)$Expected.Value
  price_down[k,1]=result
}

gamma_greek_x0=(mean(price_up[,1])+mean(price_down[,1])-2*mean(price_stable[,1]))/(delta_price^2)


Francesco Bianchi
Matricola 4506994
francesco.bianchi04 at icatt.it<mailto:francesco.bianchi04 at icatt.it>




----------------------------------------DISCLAIMER----------------------------------------

Il presente indirizzo di posta elettronica è attribuito automaticamente agli studenti e laureati dell’Università Cattolica del Sacro Cuore. L’utilizzo da parte dell’utente è personale e l’Università Cattolica del Sacro Cuore non è responsabile del contenuto delle comunicazioni.

This email address was assigned automatically as part of the registration process for students and graduates of Università Cattolica del Sacro Cuore. The use of the address is personal and Università Cattolica del Sacro Cuore accepts no liability for the content of any emails.

	[[alternative HTML version deleted]]