LongstaffpSchwartzp-mai

I asked ChatGPT to draft a code to simulate a Brownian Motion (BM), the idea is to study and write a dissertation on american option pricing, using the LS article as motivation and source of examples.

# Function to simulate Geometric Brownian Motion (GBM)

simulate_gbm <- function(S0, mu, sigma, T, n) {
  # S0: Initial stock price
  # mu: Expected return (drift)
  # sigma: Volatility
  # T: Time horizon (in years)
  # n: Number of time steps

  dt <- T / n  # Time step size
  t <- seq(0, T, by = dt) # Time vector
  W <- c(0, cumsum(sqrt(dt) * rnorm(n))) # Brownian motion (Wiener process)
  S <- S0 * exp((mu - 0.5 * sigma^2) * t + sigma * W) # GBM formula
  return(data.frame(time = t, price = S))
}

# Example usage:

S0 <- 100    # Initial stock price
mu <- 0.1     # Expected return (10% per year)
sigma <- 0.2  # Volatility (20% per year)
T <- 1        # Time horizon (1 year)
n <- 252      # Number of time steps (e.g., daily for a year)

# Simulate GBM
gbm_simulation <- simulate_gbm(S0, mu, sigma, T, n)

# Plot the simulated stock prices
plot(gbm_simulation$time, gbm_simulation$price,
     type = "l",
     xlab = "Time (years)",
     ylab = "Stock Price",
     main = "Simulated Geometric Brownian Motion")

#to plot multiple simulations

num_sims = 10; #number of simulations
all_sims = list() #list to store the simulation results.

for(i in 1:num_sims){
  all_sims[[i]] = simulate_gbm(S0, mu, sigma, T, n)
}

plot(all_sims[[1]]$time, all_sims[[1]]$price, type = "l", xlab = "Time (years)", ylab = "Stock Price", main = "Multiple GBM Simulations", ylim = range(sapply(all_sims, function(x) x$price))) #set ylim to include all simulations.

for(i in 2:num_sims){
  lines(all_sims[[i]]$time, all_sims[[i]]$price, type = "l")
}

now i tried to replicate the outcomes in Table 1 with a simulation, assuming assets are european, ie they can exercised only at maturity \(T\):

res <- replicate(10000,simulate_gbm(40,0.06,0.2,1,50)[51,2])
mean(pmax(40-res,0)*exp(-0.06))
[1] 2.093069

the correct (european) results is 2.066, as can be seen in the 9th line of Table 1. many other results can be re-obtained. for instance, the last line would approximately be using only 10000 simulations

res <- replicate(10000,simulate_gbm(44,0.06,0.4,2,50)[51,2])
mean(pmax(40-res,0)*exp(-2*0.06))
[1] 5.246892

at this stage we/you could:

but the more interesting part would be replicate Figure 1 and find the exercise boundary. this would be a bit like finding the red line below

plot(all_sims[[1]]$time, all_sims[[1]]$price, type = "l", xlab = "Time (years)", ylab = "Stock Price", main = "Multiple GBM Simulations", ylim = range(sapply(all_sims, function(x) x$price))) #set ylim to include all simulations.

for(i in 2:num_sims){
  lines(all_sims[[i]]$time, all_sims[[i]]$price, type = "l")
}

f <- function(time) 85+15*exp(-time/20)
curve(f,add=TRUE,col=2,lwd=2)

i just tried with a parametric exercise boundary such as \(a+b\exp(-dt)\), for some \(a,b,d\)

gbm_simulation <- simulate_gbm(40,0.06, 0.2, 1, 50)
# the previous one is a simulaton and we have to check if there is any early exercise
f <- function(time) 34+6*exp(-3*time)
times <- seq(0,1,len=51)
which(gbm_simulation$price-f(times)<0)
 [1]  2  3  4  5  9 10 11 12 13 14 15 16 17 18 19 22 23 24
exercise <- function(time,a=34,b=6,d=3,offset=1) a+b*exp(-d*time)-offset
discounted_payoff <- function(x,...) {
  negative <- which(x$price-exercise(x$time,...)<0)[1]
  if(is.na(negative)) negative <- 51
  dp <- max(0,40-x$price[negative])*exp(-times[negative]*0.06)
  # plot(x,t="b",main=negative);grid(col=1);curve(exercise(x,...),add=TRUE,col=2)
  dp
}
x <- simulate_gbm(40,0.06, 0.2, 1, 50)
cat("pippo\n")
pippo
set.seed(123)
res <- replicate(10000,discounted_payoff(simulate_gbm(40,0.06, 0.2, 1, 50),offset=2.5))
mean(res)
[1] 2.183727

i can create a function, optimize and plot a graph

profit <- function(offset, d,myseed=123,numsim=1000){
  set.seed(myseed)
  res <- replicate(numsim,discounted_payoff(simulate_gbm(40,0.06, 0.2, 1, 50),offset=offset,d=d))
  mean(res)
}
profitb <- function(x,...) profit(x[1],x[2],...)
profit(2.5,3)
[1] 2.125034
profitb(c(2.5,3))
[1] 2.125034
optim(c(2,3),profitb,control=list(fnscale=-1))
$par
[1] 2.155420 2.403076

$value
[1] 2.211575

$counts
function gradient 
      37       NA 

$convergence
[1] 0

$message
NULL

now the previous exercise region can be plotted against simulations, reusing and copy-pasting previous code:

num_sims = 20; #number of simulations
all_sims = list() #list to store the simulation results.

for(i in 1:num_sims){
  all_sims[[i]] = simulate_gbm(40, 0.06, 0.2, 1, 50)
}

plot(all_sims[[1]]$time, all_sims[[1]]$price, type = "l", xlab = "Time (years)", ylab = "Stock Price", main = "Multiple GBM Simulations", ylim = range(sapply(all_sims, function(x) x$price))) #set ylim to include all simulations.

for(i in 2:num_sims){
  lines(all_sims[[i]]$time, all_sims[[i]]$price, type = "l")
}
grid(col=1)
#parameters offset and d are taken from the previous optim
curve(exercise(x,offset=2.16,d=2.40),add=TRUE,col=2,lwd=2)

replotting Figure 1 in LS, i guess it’s exercise / strike (but take care: 50 vs 52 weeks, must be checked/fixed)

x <- seq(0,1,len=51)
plot(x*52,sapply(x,exercise,offset=2.16,d=2.40)/40,t="b",pch=18,ylim=c(0.7,1.05));grid(col=1)

finally, what is the price of the american option? (for the parameters that were used to estimate the optimal exercise)

during the optim, with those given 1000 paths for a given seed it was 2.21. now a run a “serious” simulation

res <- replicate(100000,discounted_payoff(simulate_gbm(40,0.06, 0.2, 1, 50),offset=2.16,d=2.40))
mean(res)
[1] 2.186576
cat("--- *** ---\n")
--- *** ---
res <- replicate(100000,discounted_payoff(simulate_gbm(40,0.06, 0.2, 1, 50),offset=3.13,d=2.60))
mean(res)
[1] 2.208117

the correct american result is 2.314, taken from Table 1. overall, our results are perfectly in line with LS and we have reasons to be happy!

Quarto enables you to weave together content and executable code into a finished document. To learn more about Quarto see https://quarto.org.

Running Code

When you click the Render button a document will be generated that includes both content and the output of embedded code. You can embed code like this:

1 + 1
[1] 2

You can add options to executable code like this

[1] 4

The echo: false option disables the printing of code (only output is displayed).