R Labs 5

Brownian motion, binomial trees and Monte Carlo simulations.

R Example 5.1 (Brownian motion): R commands to create and plot an approximate sample path of an arithmetic Brownian motion for given α and σ, over the time interval [0,T] and with n points.

alpha=0; sigma=1; T=1; n=2^(12); X0=0.1;
#############Generate 1 trajectory


R Example 5.2 (Geometric Brownian motion): For a given stock with expected rate of return μ and volatility σ, and initial price P0 and a time horizon T, simulate in R nt many trajectories of the price Pt from time  t=0  up until  t=T  through n many time periods, each of length  Δt = T/n, assuming the geometric Brownian motion model.
This is done with the function  GBM in the package  sde.

mu=0.16; sigma=0.2; P0=40; T = 1/12 ##1 month
nt=50; n=2^(8)
#############Generate nt trajectories
dt=T/n; 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,],t='l',ylim=c(ymin, ymax), col=1,
    ylab="Price P(t)",xlab="time t")
for(i in 2:nt){lines(t,X[i,], t='l',ylim=c(ymin, ymax),col=i)}


R Example 5.6 (Binomial tree option valuation): To compute the price of an American put option on a stock with current value of 50, exercise price 50, time to maturity 5 months, annualized rate of interest r is 10%, annualized volatility σ of the stock is of 40%, the annualized cost-of-carry rate b in this case equals the rate of interest r, and the number of time steps n = 5. Use methods CRRBinomialTreeOptionBinomialTreeOptionBinomialTreePlot in fOptions package from Rmetrics.


CRRBinomialTreeOption(TypeFlag = "pa", S = 50, X = 50,
     Time = 5/12, r = 0.1, b = 0.1, sigma = 0.4, n = 5)
## output : Option Price: 4.488459

##to plot an option tree (for put american)
CRRTree = BinomialTreeOption(TypeFlag = "pa", S = 50, X = 50,
      Time = 0.4167, r = 0.1, b = 0.1, sigma = 0.4, n = 5)
BinomialTreePlot(CRRTree, dy = 1, cex = 0.8, ylim = c(-6, 7),
    xlab = "n", ylab = "Option Value")
title(main = "Option Tree")

##re-run with TypeFlag = "ca" to get price for corresponding American call


R Example 5.7 (Monte Carlo option valuation):
To do a Monte Carlo simulation of arithmetic Asian option using Brownian paths with pseudo random numbers. The following step-by-step instructions for doing MC is based on the examples from the on-line manual of method MonteCarloOptions in fOptions (Rmetrics).


## MC simulation
## Step 1: function to generate the option's innovations.
## Use normal pseudo random numbers
pseudoInnovations = function(mcSteps, pathLength,init){
  #Create normal pseudo random numbers
  innovations = rnorm.pseudo(mcSteps, pathLength,init)
  #Return value
       innovations }

## Step 2: function to generate the option's price paths.
wienerPath = function(eps) {
  # Generate the Paths:
  path = (b-sigma*sigma/2)*delta.t + sigma*sqrt(delta.t)*eps
  # Return Value:
       path }

## Step 3: function for the payoff for an Asian Call or Put:
arithmeticAsianPayoff = function(path) {
  # Compute the Call/Put Payoff Value:
  SM = mean(S*exp(cumsum(path)))
  if (TypeFlag == "c") payoff = exp(-r*Time)*max(SM-X, 0)
  if (TypeFlag == "p") payoff = exp(-r*Time)*max(0, X-SM)
  # Return Value:
       payoff }

## Final Step: Set global parameters for the  option:
TypeFlag <<- "c"; S <Time <
# Do the Asian Simulation with pseudo random numbers:
mc = MonteCarloOption(delta.t=1/360, pathLength=30, mcSteps=5000,
  mcLoops = 50, init = TRUE, innovations.gen = pseudoInnovations,
  path.gen = wienerPath, payoff.calc = arithmeticAsianPayoff,
  antithetic = TRUE, standardization = FALSE, trace = TRUE)
# Plot the MC Iteration Path:
par(mfrow = c(1, 1))
mcPrice = cumsum(mc)/(1:length(mc))
plot(mcPrice, type = "l", main = "Arithmetic Asian Option",
     xlab = "Monte Carlo Loops", ylab = "Option Price")

To compare the performance of the MC method against other analytical methods, consider for example the Turnbull and Wakeman's approximation method.

    X=100, Time=1/12, time=1/12, tau=0, r=0.1, b=0.1, sigma=0.4)