incorporating_measurement.R

less_uncertainty_smaller_step <- 
  .25 * (1:1000 %>% map(~ run_simulation(initial = 1 - 0.1)) %>% unlist %>% mean) + 
    .5 * (1:1000 %>% map(~ run_simulation(initial = 1)) %>% unlist %>% mean) + 
    .25 * (1:1000 %>% map(~ run_simulation(initial = 1 + 0.1)) %>% unlist %>% mean) 
  
  

greater_uncertainty_bigger_step <- 
  .05 * (1:1000 %>% map(~ run_simulation(initial = 1 - 3, step_size = 0.5)) %>% unlist %>% mean) + 
  .1 * (1:1000 %>% map(~ run_simulation(initial = 1 - 2, step_size = 0.5)) %>% unlist %>% mean) + 
  .2 * (1:1000 %>% map(~ run_simulation(initial = 1 - 1, step_size = 0.5)) %>% unlist %>% mean) + 
  .3 * (1:1000 %>% map(~ run_simulation(initial = 1, step_size = 0.5)) %>% unlist %>% mean) + 
  .2 * (1:1000 %>% map(~ run_simulation(initial = 1 + 1, step_size = 0.5)) %>% unlist %>% mean) + 
  .1 * (1:1000 %>% map(~ run_simulation(initial = 1 + 2, step_size = 0.5)) %>% unlist %>% mean) + 
  .05 * (1:1000 %>% map(~ run_simulation(initial = 1 + 3, step_size = 0.5)) %>% unlist %>% mean)

option_pricing_under_uncertainty.R

library(e1071)
library(tidyverse)
run_simulation <- function(
  step_size = 0.1, 
  steps = 100, 
  initial = 1, 
  min_val = -5, 
  max_val = 5) {
  
  add_step <- function(val) {
    
    candidate_val <- val + sample(c(0, step_size, -step_size), 1)
    if(candidate_val > max_val) return(max_val)
    if(candidate_val < min_val) return(min_val)
    return(candidate_val)
    
  }
  
  current_val <- initial
  for(i in 1:steps) {
    current_val <- add_step(current_val)
  } 
  
  sig <- sigmoid(current_val)
  
  result <- runif(1) < sig
  
  return(result)
  
}

sigmoid(1)
1:1000 %>% map(~ run_simulation()) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation(step_size = 0.3)) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation(step_size = 0.5)) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation(step_size = 1)) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation(step_size = 5)) %>% unlist %>% mean



run_simulation_without_constraints <- function(
  step_size = 0.1, 
  steps = 100, 
  initial = 1) { 
  # min_val = -5, 
  # max_val = 5
  
  add_step <- function(val) {
    
    candidate_val <- val + sample(c(0, step_size, -step_size), 1)
    # if(candidate_val > max_val) return(max_val)
    # if(candidate_val < min_val) return(min_val)
    return(candidate_val)
    
  }
  
  current_val <- initial
  for(i in 1:steps) {
    current_val <- add_step(current_val)
  } 
  
  sig <- sigmoid(current_val)
  
  result <- runif(1) < sig
  
  return(result)
  
}


sigmoid(1)
1:1000 %>% map(~ run_simulation_without_constraints()) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation_without_constraints(step_size = 0.3)) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation_without_constraints(step_size = 0.5)) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation_without_constraints(step_size = 1)) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation_without_constraints(step_size = 5)) %>% unlist %>% mean

run_simulation_on_sigmoid <- function(
  step_size = 0.1, 
  steps = 100, 
  initial = 1, 
  min_val = -5, 
  max_val = 5) {
  
  add_step <- function(val) {
    
    candidate_val <- val + sample(c(0, step_size, -step_size), 1)
    if(candidate_val > max_val) return(max_val)
    if(candidate_val < min_val) return(min_val)
    return(candidate_val)
    
  }
  
  current_val <- initial
  for(i in 1:steps) {
    current_val <- add_step(current_val)
  } 
  
  sig <- sigmoid(current_val)
  
  
  
  return(sig)
  
}


sigmoid(1)
1:1000 %>% map(~ run_simulation_on_sigmoid()) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation_on_sigmoid(step_size = 0.3)) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation_on_sigmoid(step_size = 0.5)) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation_on_sigmoid(step_size = 1)) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation_on_sigmoid(step_size = 5)) %>% unlist %>% mean


run_simulation_on_sigmoid_without_constraints <- function(
  step_size = 0.1, 
  steps = 100, 
  initial = 1) { 
  # min_val = -5, 
  # max_val = 5) {
  
  add_step <- function(val) {
    
    candidate_val <- val + sample(c(0, step_size, -step_size), 1)
    # if(candidate_val > max_val) return(max_val)
    # if(candidate_val < min_val) return(min_val)
    return(candidate_val)
    
  }
  
  current_val <- initial
  for(i in 1:steps) {
    current_val <- add_step(current_val)
  } 
  
  sig <- sigmoid(current_val)
  
  
  
  return(sig)
  
}


sigmoid(1)
1:1000 %>% map(~ run_simulation_on_sigmoid_without_constraints()) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation_on_sigmoid_without_constraints(step_size = 0.3)) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation_on_sigmoid_without_constraints(step_size = 0.5)) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation_on_sigmoid_without_constraints(step_size = 1)) %>% unlist %>% mean
1:1000 %>% map(~ run_simulation_on_sigmoid_without_constraints(step_size = 5)) %>% unlist %>% mean

simulated_price_paths.R

run_simulation_of_sim_paths_on_sigmoid_with_constraints <- function(
  step_size = 0.1, 
  steps = 100, 
  initial = 1,
  sims_per_step = 100, 
  max_val = 5, 
  min_val = -5
  ) { 
  
  
  add_step <- function(val) {
    
    candidate_val <- val + sample(c(0, step_size, -step_size), 1)
    if(candidate_val > max_val) return(max_val)
    if(candidate_val < min_val) return(min_val)
    return(candidate_val)
    
  }
  
  current_val <- initial
  current_prob <- c()
  price_list <- c()
  current_price <- c()
  
  for(i in 0:steps) {
    
    current_val <- add_step(current_val)
    remaining_steps <- steps-i
    
    if(remaining_steps > 0) {
      
      simulated_outcomes <- c()
      for(j in 1:sims_per_step) {
        
        sim_val <- current_val
        for(k in 1:remaining_steps){
          sim_val <- add_step(sim_val)
        }
        simulated_outcomes <- c(simulated_outcomes, sim_val)
      }
      current_price <- c(current_price, mean(sapply(simulated_outcomes, sigmoid)))
      current_prob <- c(current_prob, sigmoid(current_val))
    }
    
    
  } 
  
  
  
  return(data.frame(price = current_price, prob = current_prob))
  
}

run_simulation_of_sim_paths_on_sigmoid_without_constraints <- function(
  step_size = 0.1, 
  steps = 100, 
  initial = 1,
  sims_per_step = 100
) { 
  
  
  add_step <- function(val) {
    
    candidate_val <- val + sample(c(0, step_size, -step_size), 1)
    # if(candidate_val > max_val) return(max_val)
    # if(candidate_val < min_val) return(min_val)
    return(candidate_val)
    
  }
  
  current_val <- initial
  current_prob <- c()
  price_list <- c()
  current_price <- c()
  
  for(i in 0:steps) {
    
    current_val <- add_step(current_val)
    remaining_steps <- steps-i
    
    if(remaining_steps > 0) {
      
      simulated_outcomes <- c()
      for(j in 1:sims_per_step) {
        
        sim_val <- current_val
        for(k in 1:remaining_steps){
          sim_val <- add_step(sim_val)
        }
        simulated_outcomes <- c(simulated_outcomes, sim_val)
      }
      current_price <- c(current_price, mean(sapply(simulated_outcomes, sigmoid)))
      current_prob <- c(current_prob, sigmoid(current_val))
    }
    
    
  } 
  
  
  
  return(data.frame(price = current_price, prob = current_prob))
  
}


.t1 <- run_simulation_of_sim_paths_on_sigmoid_with_constraints(initial = 0.5, step_size = 0.1, steps = 300, sims_per_step = 100) 
.t2 <- run_simulation_of_sim_paths_on_sigmoid_with_constraints(initial = 0.5, step_size = 0.4, steps = 300, sims_per_step = 100)


.t1 %>% 
  ggplot() + 
  aes(x = 1:nrow(.t1)) + 
  geom_line(aes(y = price), color = "red") + 
  geom_line(aes(y = prob), color = "blue") 

.t2 %>% 
  ggplot() + 
  aes(x = 1:nrow(.t2)) + 
  geom_line(aes(y = price), color = "red") + 
  geom_line(aes(y = prob), color = "blue") 


library(foreach)
library(doParallel)
registerDoParallel(4)

non_vol_runs <- foreach(i = 1:100, .combine = bind_rows) %dopar% {
  sim <- run_simulation_of_sim_paths_on_sigmoid_with_constraints(initial = 0, step_size = 0.1, steps = 600, sims_per_step = 100) 
  sim %>% mutate(
    run = i,
    type = "non_vol"
  )
}

vol_runs <- foreach(i = 1:100, .combine = bind_rows) %dopar% {
  sim <- run_simulation_of_sim_paths_on_sigmoid_with_constraints(initial = 0, step_size = 0.4, steps = 600, sims_per_step = 100) 
  sim %>% mutate(
    run = i,
    type = "vol"
  )
}

vol_runs %>% 
  filter(run == 2) %>% 
  ggplot() + 
  aes(x = 1:300) + 
  geom_line(aes(y = price), color = "red") + 
  geom_line(aes(y = prob), color = "blue") 

saveRDS(non_vol_runs, "~/Downloads/non_vol_runs.RDS")
saveRDS(vol_runs, "~/Downloads/vol_runs.RDS")

bind_rows(non_vol_runs, vol_runs) %>% 
  group_by(type, run) %>% 
  mutate(idx = 1:n()) %>% 
  gather(key = series, value = value, -idx, -run, -type) %>% 
  # filter(run %in% c(1, 2, 3)) %>% 
  mutate(grp = str_c(run, "_", series)) %>% 
  filter(series == "price") %>% 
  ggplot() + 
  aes(x = idx) + 
  aes(y = value) +
  # aes(color = series) +
  geom_line(aes(group = grp), alpha = 0.5) +
  labs(y = "Option price", x = "Time") + 
  # geom_density_2d(aes(group = series)) + 
  facet_wrap(~ type) + 
  theme_bw()

bind_rows(non_vol_runs, vol_runs) %>% 
  group_by(type, run) %>% 
  mutate(idx = 1:n()) %>% 
  gather(key = series, value = value, -idx, -run, -type) %>% 
  group_by(idx, series, type) %>% 
  summarize(
    # quant_99 = quantile(value, .99),
    quant_95 = quantile(value, .90),
    # quant_80 = quantile(value, .8),
    # quant_90 = quantile(value, .9),
    # quant_98 = quantile(value, .98),
    quant_05 = quantile(value, .05),
    # quant_02 = quantile(value, 0.02)
    # quant_01 = quantile(value, 0.01)
  ) %>% 
  ungroup %>% 
  gather(key = qntle, value = value, -series, -idx, -type) %>% 
  mutate(grp = str_c(qntle, "_", series, "_", type)) %>% 
  ggplot() + 
  aes(x = idx) + 
  aes(y = value) + 
  aes(group = grp) + 
  aes(color = type) +
  aes(linetype = series) +
  geom_smooth(se = F, lwd = 0.7) + 
  scale_color_brewer("", type = "qual", labels = c("Step size = 0.1", "Step size = 0.4")) + 
  scale_linetype_discrete("", labels = c("Price of option", "Underlying probability")) + 
  guides(lty = guide_legend(override.aes = list(col = 'darkgrey'))) + 
  theme_bw() + 
  theme(legend.position = c(0.6, 0.54), 
        legend.background = element_rect(fill = "transparent"),
        legend.box = "horizontal", 
        legend.key = element_rect(fill = "transparent")) + 
  labs(x = "Time", y = "Value", title = "Trading bounds for different volatilities", 
       caption = "5% and 95% quantile, respectively\n100 runs for each step size\nPrice computed from 100 simulations at each step\nUnderlying bound in [-5,5]\nOutcome binomial draw from sigmoid transformation") +  
  # facet_wrap(~ type) + 
  NULL