Estimate barrier breach probabilities via simulation, coded into an RMarkdown document

barrier-breach-simulation.Rmd

---
title: "Barrier Breach Simulation"
author: "Paul Teetor"
date: "June 2022"
output: html_document
params:
  annualVol: 0.20
  annualDrift: 0.0
  term: 30   # In days
  startPrice: 100
  lowerBound: 90.9
  upperBound: 110
  ndraws: 5000
---

```{r setup, include=FALSE}
suppressPackageStartupMessages({
  library(tidyverse)
})

knitr::opts_chunk$set(echo = FALSE)
```

```{r parameters}
annualVol = as.double(params$annualVol)
annualDrift = as.double(params$annualDrift)
term = as.integer(params$term)   # days

startPrice = as.double(params$startPrice)
upperBound = as.double(params$upperBound)
lowerBound = as.double(params$lowerBound)
ndraws = as.integer(params$ndraws)

stepVol = annualVol / sqrt(252)
stepDrift = annualDrift / 252


```

```{r functions, include = FALSE}
onePath = function() {
  r <- rnorm(term, mean = stepDrift, sd = stepVol)
  startPrice * cumprod(exp(c(0, r)))
}

smartMax = function(x, y) ifelse(is.na(x) & is.na(y), NA,
                                 pmax(x, y, na.rm = TRUE))

smartMin = function(x, y) ifelse(is.na(x) & is.na(y), NA,
                                 pmin(x, y, na.rm = TRUE))

oneDraw = function(i) {
  path <- onePath()
  
  xlower <- which(path <= lowerBound)[1]   # NA if never crosses
  xupper <- which(path >= upperBound)[1]   # NA if never crosses
  
  tibble(
    Draw = i,
    CrossLower = !is.na(xlower),
    CrossLowerTime = xlower,
    CrossUpper = !is.na(xupper),
    CrossUpperTime = xupper,
    CrossEither = CrossLower | CrossUpper,
    CrossEitherTime = smartMin(CrossLowerTime, CrossUpperTime)
  )
}
```

```{r simulation}
draws <- ((1:ndraws)
          |> map_dfr(oneDraw) )

simSummary <-
  with(draws,
       tibble(ProbCrossLower = mean(CrossLower),
              MeanTimeToCrossLower = mean(CrossLowerTime, na.rm = TRUE),
              ProbCrossUpper = mean(CrossUpper),
              MeanTimeToCrossUpper = mean(CrossUpperTime, na.rm = TRUE),
              ProbCrossEither = mean(CrossEither),
              MeanTimeToCrossEither = mean(CrossEitherTime, na.rm = TRUE) )
  )
```

* Annualized volatility: `r scales::percent(annualVol, accuracy = 0.1)`
* Annualized drift: `r scales::percent(annualDrift, accuracy = 0.1)`
* Term: `r term` days
* Starting price: `r startPrice`
* Lower bound: `r lowerBound`
* Upper bound: `r upperBound`
* Number of simulated paths: `r scales::comma(ndraws)`

<br/>

```{r show-prob-summary}
(simSummary
 |> select(starts_with("Prob"))
 |> gt::gt()
 |> gt::tab_header(title = "Probability of breach")
 |> gt::fmt_percent(starts_with("Prob"), decimals = 1)
 |> gt::cols_label(ProbCrossLower = "Lower breach",
                   ProbCrossUpper = "Upper breach",                        
                   ProbCrossEither = "Either" ))
```

<br/>

```{r show-time-summary}
(simSummary
 |> select(starts_with("MeanTime"))
 |> gt::gt()
 |> gt::tab_header(title = "Mean time to breach")
 |> gt::fmt_number(starts_with("MeanTime"), decimals = 1)
 |> gt::cols_label(MeanTimeToCrossLower = "Lower breach",
                   MeanTimeToCrossUpper = "Upper breach",
                   MeanTimeToCrossEither = "Either" ))
```

<br/>
The *mean time to breach* estimates are conditional upon the event occuring
and do not include right-censored data.
<br/>
<br/>

```{r plot-distributions}
caption <- paste0("Annual vol = ", annualVol, ", annual drift = ", annualDrift)
subtitle <- paste0("lower bound = ", lowerBound, " , upper bound = ", upperBound)

p <- (draws
      |> dplyr::filter(!is.na(CrossLowerTime) | !is.na(CrossUpperTime))
      |> with(rbind(tibble(Bound = "Lower", Time = CrossLowerTime),
                    tibble(Bound = "Upper", Time = CrossUpperTime) ))
      |> ggplot(aes(x = Time, color = Bound))
      + geom_density()
      + labs(title = "Distribution of time to breach",
             subtitle = subtitle,
             caption = caption ) )

suppressWarnings(print(p))
```