3

Simulate in a graph 50 sample paths of a stock price $80 over 90 days modeled as geometric Brownian motion with drift parameter 0.1 and volatility 0.5. Show in a graph this process on the vertical axis Price option and time on the horizontal axis. Find the probability that in 90 days the option price of will rise to at least $100.

library(sde)
mu<-0.1
sigma<-0.5
P0<-80 #initial price
T <-90/360 #time in years
nt=10000 #number of trajectories within each simulation
n=100 #time periods
dt<-T/n #length of time periods
t <- seq(0,T,by=dt)
X=matrix(rep(0,length(t)*nt),nrow = nt)
for(i in 1:nt) {
  X[i,]=GBM(x=P0, r=mu, sigma=sigma, T=T, N=n)
}
ymax=max(X); ymin=min(X) #bounds for simulated prices
plot(t,X[1,],type="l", ylim=c(ymin,ymax), col=1, xlab="Time", ylab="Price Y(t)")
for(i in 2:nt){
  lines(t,X[i,], type='l', ylim=c(ymin, ymax), col=i)
}
Prob<-sum(nt>=100)/nt
Prob
desertnaut
  • 57,590
  • 26
  • 140
  • 166
Andrea Garcia
  • 127
  • 2
  • 8
  • 3
    So you simulated 50 stock paths. How many of them ended up above the target price of $100? – rajah9 May 09 '20 at 12:24
  • 2
    50 sample paths is not enough to get a good estimate of the probability. You could perhaps graph 50, but then estimate the probability with e.g. 10,000 paths. (Better yet both give the estimate based on the 50 shown, and also report the results of a larger simulation). – John Coleman May 09 '20 at 12:26
  • OK. I changed the number of time periods to 10000 but how could do it using Monte carlo simulation? – Andrea Garcia May 09 '20 at 12:42
  • Why change the number of time periods? It should be the number of paths through that time period which should be changed. In any event, once you do the simulation, what is the problem? Just count the paths that end up above the target prices and divide by the number of paths. Monte Carlo simulation is fundamentally a very naive algorithm. If you want to estimate the probability of rolling a seven in a pair of dice, just roll it 100 times, count the sevens, and divide by 100. That is all that is happening here, except you have stock price paths rather than dice. – John Coleman May 09 '20 at 12:59
  • 10,000 sample paths looks like a bunch of trajectories in the graph and I only need 50. Also I tried to estimate the probability of at least $100 and something has to be wrong. Thanks a lot! – Andrea Garcia May 09 '20 at 15:01
  • 1
    The simulation itself and what you do with that simulation are different things. You could graph it, estimate a probability with it, etc. Personally, I would write a function which simulates a single path, debug that function, then `replicate` it any number of times. As far as your code goes, `sum(nt>=100)/nt` doesn't make sense. `nt` is just a single number so what is the point of `sum(nt>=100)`? You are not summing the right thing. – John Coleman May 09 '20 at 15:19
  • Could you please explain it tp me in detail? An example about how to simulate a single path and then replicated it. What should I sum in order to get the probability? – Andrea Garcia May 09 '20 at 17:55

1 Answers1

1

The answer depends on how you interpret your parameters t, n and T.

I make a few assumptions: since

T <- 90/360

i suppose, this means 90 days of a year (approx. 360, which is common in finance for a year). Your definition of t

n <- 100
dt <- T/n
t <- seq(0, T, by=dt)

gives the timesteps for your simulation, therefore your day 90 is simply given by max(t) = 0.25 = T with index 101, the last element of t.

X contains 50 paths of your stochastic process, indexed by X[i,], at timesteps j given by X[,j]. So if you want to know the values of your 50 simulations at day 90 just look at X[,101].

So you want to know how many of your paths exceed 100 at day 90. Just count them by

success <- sum(X[,101] >= 100)

If you want to calculate a empirical probability just divide them by the number of your paths. Therefore

emp_prob <- success/nt
Martin Gal
  • 16,640
  • 5
  • 21
  • 39
  • Thanks so much @Martin Gal. How I could estimate the maximum of an option price will rise to at least 100? – Andrea Garcia May 11 '20 at 13:47
  • Try something like `success <- sum(sapply(1:n, function(i) max(X[i,]) >=100))` and then calculate the empirical probability like shown above. – Martin Gal May 11 '20 at 14:00