# R/JdmbsNewJump.R In Jdmbs: Monte Carlo Option Pricing Algorithms for Jump Diffusion Models with Correlational Companies

#### Documented in jdm_new_bs

```#' A Monte Carlo Option Pricing Algorithm for Jump Diffusion Model with Correlation Companies
#' @param  companies_data : a matrix of a correlation coefficient of companies
#' @param  companies : an integer of a company number in order to simulate.
#' @param  simulation.length : an integer of a time duration of simulation.
#' @param  monte_carlo : an integer of an iteration number for monte carlo.
#' @param  start_price : a vector of company's initial stock prices.
#' @param  mu : a vector of drift parameters of geometric Brownian motion.
#' @param  sigma : a vector of volatility parameters of geometric Brownian motion.
#' @param  event_times : an integer of how many times jump in unit time.
#' @param  jump : a vector of jump parameter.
#' @param  K : a vector of option strike prices.
#' @param  color : a vector of colors in plot.
#' @return option prices : a list of (call_price, put_price)
#' @examples
#' price <- jdm_new_bs(matrix(c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),nrow=3, ncol=3),
#'                     3, simulation.length=100,monte_carlo=80, c(1000,500,500),
#'                     c(0.002, 0.012, 0.005),c(0.05,0.05,0.06), 3,c(0.1,0.1,0.1),
#'                     c(1500,1000,700),c("red","blue","green")
#'                    )
#' @export
jdm_new_bs<- function(companies_data, companies, simulation.length=180, monte_carlo=1000, start_price=start_price, mu=mu,sigma=sigma, event_times=event_times, jump=jump, K=K, color=color) {
if(is.matrix(companies_data) == FALSE){
print("Error: Input companies_data type is not matrix.")
return(FALSE)
}
if(is.numeric(companies) == FALSE){
print("Error: Input companies type is not numeric.")
return(FALSE)
}
if(companies <= 0){
print("Error: Input companies is less than or equal to zero.")
return(FALSE)
}
if(is.numeric(simulation.length) == FALSE){
print("Error: Input simulation.length type is not integer.")
return(FALSE)
}
if(simulation.length <= 0){
print("Error: Input simulation.length is less than or equal to zero.")
return(FALSE)
}
if(is.numeric(monte_carlo) == FALSE){
print("Error: Input monte_carlo type is not integer.")
return(FALSE)
}
if(monte_carlo <= 0){
print("Error: Input monte_carlo is less than or equal to zero.")
return(FALSE)
}
price_sim <-array(0, c(companies, monte_carlo));
price <- array(0, c(monte_carlo, companies, simulation.length))
price_up_limit_graph <- 0;
price_down_limit_graph <- 0;
dW <- array(0, c(companies))
W <- array(0, c(companies))
J <- array(0, c(companies))
jump_time <- array(0, c(monte_carlo, simulation.length))
jump_time_count <- array(0, c(monte_carlo))
jump_count <- 1
jump_flag <- "off"

for(count in 1 : monte_carlo){
jump_count <- 1
J <- array(0, c(companies))
J_affect <- array(0, c(companies))
jump_company_id <- c(0)
W <- array(0, c(companies))
# somethings happen 3times in 100 mins. lambda = 3 / 100
lambda = event_times / simulation.length
while(jump_time[count, jump_count] < simulation.length){
# http://preshing.com/20111007/how-to-generate-random-timings-for-a-poisson-process/ from The Art of Computer Programming
jump_time[count, jump_count + 1] = jump_time[count, jump_count] - log(1 - runif(1,min=0,max=1)) / lambda
# which company occur jump. this company's jump affect other companies. company id is memorize.
jump_company_id[jump_count + 1] <- as.integer(runif(1,min=1,max=companies+1))
J[jump_count+ 1] <- runif(1,min=0,max=1) * rnorm(1, mean = 0, sd = 1)
jump_count <- jump_count + 1
jump_time_count[count] <- jump_time_count[count] + 1
}
# jump_count <- 2 because jump_time[count, 1] = 0. So start from "2"
jump_count <- 2
for(i in 2 : simulation.length){
for(j in 1 : companies){
if(i == 2){
price[count, j, 1] <- start_price[j]
}
if(jump_time[count, jump_count] <= i){
J_affect[j] <- J_affect[j] + J[jump_count] * companies_data[jump_company_id[jump_count],j]
jump_flag <- "on"
}
# W2 - W1 = N(mean=0, sd = 2 - 1)
dW[j] <- rnorm(1, mean = 0, sd = 1)
W[j] <- W[j] + dW[j]
# about Geometric Brownian Motion from Wikipedia ( https://en.wikipedia.org/wiki/Geometric_Brownian_motion )
price[count, j,i] <- start_price[j] * exp((mu[j] - (sigma[j] ^ 2) / 2) * (i - 1)  + sigma[j] * W[j] + jump[j] * J_affect[j])
if(price[count, j,i] > price_up_limit_graph){
price_up_limit_graph <- price[count, j,i]
}
if(price[count, j,i] < price_down_limit_graph){
price_down_limit_graph <- price[count, j,i]
}
}
if(jump_flag == "on"){
jump_count <- jump_count + 1
jump_flag = "off"
}
}
for (i in 1 : companies){
price_sim[i,count] <- price[count, i, simulation.length];
}
}
for(count in 1 : monte_carlo){
if(count == 1){
plot(price[count, 1, ], xlab="t", ylab="price", ylim = c(price_down_limit_graph, price_up_limit_graph), xlim=c(1, simulation.length), xaxp=c(0, simulation.length,simulation.length/10), type="l", col="black", lwd = 0.10, cex.axis = 0.6, cex.lab = 0.8)
par(new=T)
plot(x = jump_time[count,1:jump_time_count[count]], y = 1 : jump_time_count[count] - 1, type = "s", xlim = c(1, simulation.length),axes=FALSE,lwd = 0.10,xlab="",ylab="")
par(new=T)
for (i in 1 : companies){
plot(price[count, i, ], ylim=c(price_down_limit_graph, price_up_limit_graph), xlim=c(1,simulation.length), xlab="", ylab="", axes=FALSE, type="l", col=color[i], lwd = 0.20, cex.axis = 0.6)
par(new=T)
}
}
else{
plot(x = jump_time[count,1:jump_time_count[count]], y = 1 : jump_time_count[count] - 1, type = "s", xlim = c(1, simulation.length),axes=FALSE,lwd = 0.10,xlab="",ylab="")
par(new=T)
for (i in 1:companies){
plot(price[count, i, ], ylim=c(price_down_limit_graph, price_up_limit_graph), xlim=c(1,simulation.length), xlab="", ylab="", axes=FALSE, type="l", col=color[i], lwd = 0.20, cex.axis = 0.6)
par(new=T)
}
}
}
par(new=T)
title(main = "Geometric Brownian Motion", col.main="black", font.main=1, cex.main = 0.8)
#legend(1,8000, c("mu = 0.7,  sigma = 0.5"), cex=0.8, col=c("black"), pch=1, lty=1);

call_price <- array(0,c(companies));
put_price <- array(0,c(companies));
for(i in 1 : companies){
for(count in 1 : monte_carlo){
if((price_sim[i,count] - K[i]) <= 0){
call_price[i] <- call_price[i] + 0
}
else{
call_price[i] <- call_price[i] + (price_sim[i,count] - K[i])
}
}
call_price[i] <- call_price[i] / monte_carlo;
}
for(i in 1:companies){
for(count in 1 : monte_carlo){
if((K[i] - price_sim[i,count]) <= 0){
put_price[i] <- put_price[i] + 0;
}
else{
put_price[i] <- put_price[i] + (K[i] - price_sim[i,count])
}
}
put_price[i] <- put_price[i] / monte_carlo;
}
print("Call Option Price:")
print(call_price)
print("Put Option Price:")
print(put_price)
names(call_price) <- "Call Option Price"
names(put_price) <- "Put Option Price"
price <- list(call_price,put_price)
return (invisible(price))
}
```

## Try the Jdmbs package in your browser

Any scripts or data that you put into this service are public.

Jdmbs documentation built on May 2, 2018, 1:04 a.m.