EXECUTIVE SUMMARY

The below code creates functions and plots for pricing European options, American options, Asian options, Binomial models/trees, and straddles.

Note : The code uses predetermine inputs for all models.

CONTENTS

i. Pricing European Options & Binomial Models/Trees
ii. Straddle Pricing
iii.Pricing American Options
iv. Discrete Dividends
v. Pricing Asian Options by Monte Carlo
vi. Pricing American Options using the Longstaff-Schwartz LeastSquare Method

Pricing European Options & Binomial Models/Trees

Using R to write a computer code which takes as inputs:

  • The initial stock price \(S_0\)

  • The payoff function \(F(S_T)\)

  • The interest rate \(r\)

  • The length of the period \(h\)

  • The up and down factors \(u\) and \(d\)

  • The number of periods \(N\)

and which calculates the European option price as well as the composition of the replicating portfolio at every node of the tree.

stock_tree = function(S, N, u, d){
   tree = matrix(0, nrow=N+1, ncol=N+1)
   
   for(i in 1:(N+1)){
      for(j in 1:i){
         tree[i,j] = S * u^(j-1) * d^((i-1)-(j-1))
      }
   }
   tree
}

p = function(r, delta, h, u, d){
   (exp((r - delta)*h) - d)/(u - d)
}

European_option = function(S, K, N, h, r, delta, sigma, type){
   u = exp((r - delta)*h + sigma*sqrt(h))
   d = exp((r - delta)*h - sigma*sqrt(h))
   stock_tree = stock_tree(S, N, u, d)
   p = p(r, delta, h, u, d)
   option_tree = matrix(0, nrow=nrow(stock_tree), ncol=ncol(stock_tree))
   
   #Terminal nodes
   if(type == "call"){
      option_tree[nrow(option_tree),] = pmax(stock_tree[nrow(stock_tree),]-K, 0)
   }
   else if(type == "put"){
      option_tree[nrow(option_tree),] = pmax(K-stock_tree[nrow(stock_tree),], 0)
   }
   
   #Backward operations
   for(i in (nrow(option_tree)-1):1){
      for(j in 1:i){
         option_tree[i,j] = exp(-(r * h)) * (p * option_tree[i+1,j+1] + (1-p) * option_tree[i+1,j])
      }
   }
   option_tree
}

RP_stock_tree_E = function(S, K, N, h, r, delta, sigma, type){
   u = exp((r - delta)*h + sigma*sqrt(h))
   d = exp((r - delta)*h - sigma*sqrt(h))
   stock_tree = stock_tree(S, N, u, d)
   option_tree = European_option(S, K, N, h, r, delta, sigma, type)
   
   RP_stock_tree = matrix(0, nrow=nrow(stock_tree), ncol=ncol(stock_tree))
   
   for(i in (nrow(option_tree)-1):1){
      for(j in 1:i){
         RP_stock_tree[i,j] = exp(-(delta * h)) * (option_tree[i+1,j+1] - option_tree[i+1,j]) / (stock_tree[i+1,j+1] - stock_tree[i+1,j])
      }
   }
   RP_stock_tree
}

RP_bond_tree_E = function(S, K, N, h, r, delta, sigma, type){
   u = exp((r - delta)*h + sigma*sqrt(h))
   d = exp((r - delta)*h - sigma*sqrt(h))
   stock_tree = stock_tree(S, N, u, d)
   option_tree = European_option(S, K, N, h, r, delta, sigma, type)
   
   RP_bond_tree = matrix(0, nrow=nrow(stock_tree), ncol=ncol(stock_tree))
   
   for(i in (nrow(option_tree)-1):1){
      for(j in 1:i){
         RP_bond_tree[i,j] = exp(-(r * h)) * (option_tree[i+1,j] * u - option_tree[i+1,j+1] * d) / (u - d)
      }
   }
   RP_bond_tree
}

Straddle Pricing

Computing the initial value of a straddle with \(N = 4, r = 0.02, h = 0.25, u = e^{rh+0.2\sqrt{h}}, d = e^{rh-0.2\sqrt{h}}, S_0 = 100\), and strike \(K = 90\).

A straddle consists of one long call and one long put:

European_call = European_option(100, 90, 4, 0.25, 0.02, 0, 0.2, "call")
European_put = European_option(100, 90, 4, 0.25, 0.02, 0, 0.2, "put")
straddle_price = European_call[1,1] + European_put[1,1]
straddle_price
## [1] 18.48912

So the initial value of the straddle with \(N=4\) would be 18.48912.

Applying code to compute the initial value of a straddle with \(N = 40, r = 0.02, h = 0.25, u = e^{rh+0.2\sqrt{h}}, d = e^{rh-0.2\sqrt{h}}, S_0 = 100\), and strike \(K = 90\)._**

European_call = European_option(100, 90, 40, 0.25, 0.02, 0, 0.2, "call")
European_put = European_option(100, 90, 40, 0.25, 0.02, 0, 0.2, "put")
straddle_price = European_call[1,1] + European_put[1,1]
straddle_price
## [1] 48.02893

So the initial value of the straddle with \(N=40\) would be 48.02893.

Applying code to compute the initial value of a binary call option with \(N = 4, r = 0.02, h = 0.25, u = e^{rh+0.2\sqrt{h}}, d = e^{rh-0.2\sqrt{h}}, S_0 = 100\), and strike \(K = 90\)._**

Payoff of binary options is 100 if \(S_T>K\) and 0 otherwise:

binary_call = function(S, K, N, h, r, delta, sigma){
   u = exp((r - delta)*h + sigma*sqrt(h))
   d = exp((r - delta)*h - sigma*sqrt(h))
   stock_tree = stock_tree(S, N, u, d)
   p = p(r, delta, h, u, d)
   binary_call_tree = matrix(0, nrow=nrow(stock_tree), ncol=ncol(stock_tree))
   
   #Terminal nodes
   binary_call_tree[nrow(binary_call_tree),] = unlist(lapply(stock_tree[nrow(stock_tree),], function(x) if(x > K) {100} else{0}))
   
   #Backward operations
   for(i in (nrow(binary_call_tree)-1):1){
      for(j in 1:i){
         binary_call_tree[i,j] = exp(-(r * h)) * (p * binary_call_tree[i+1,j+1] + (1-p) * binary_call_tree[i+1,j])
      }
   }
   binary_call_tree
}

binary_call = binary_call(100, 90, 4, 0.25, 0.02, 0, 0.2)
binary_call
##          [,1]     [,2]      [,3]      [,4] [,5]
## [1,] 63.62740  0.00000   0.00000   0.00000    0
## [2,] 45.56757 84.25801   0.00000   0.00000    0
## [3,] 22.33996 71.71890  99.00498   0.00000    0
## [4,]  0.00000 47.26516  99.50125  99.50125    0
## [5,]  0.00000  0.00000 100.00000 100.00000  100

So the initial value of the binary call would be 63.62740.

Pricing American Options

Computer code which takes as inputs:

  • The initial stock price \(S_0\)

  • The payoff function \(g(S_T)\)

  • The interest rate \(r\)

  • The length of the period \(h\)

  • The up and down factors \(u\) and \(d\)

  • The number of periods \(N\)

and which calculates the American option price as well as the composition of the replicating portfolio at every node of the tree and also determines the optimal exercise dates.

The initial value of an American put and an American call with strike \(K = 10\) in a binomial model with \(r = 0.01, u = e^{rh+0.15\sqrt{h}}, d = e^{rh-0.15\sqrt{h}}, h = 1/365, S_0 = 10\), and \(N =250\) periods._**

American_option = function(S, K, N, h, r, delta, sigma, type){
   u = exp((r - delta)*h + sigma*sqrt(h))
   d = exp((r - delta)*h - sigma*sqrt(h))
   stock_tree = stock_tree(S, N, u, d)
   p = p(r, delta, h, u, d)
   option_tree = matrix(0, nrow=nrow(stock_tree), ncol=ncol(stock_tree))
   
   if(type == "call"){
      #Terminal nodes
      option_tree[nrow(option_tree),] = pmax(stock_tree[nrow(stock_tree),]-K, 0)
      
      #Backward operations
      for(i in (nrow(option_tree)-1):1){
         for(j in 1:i){
            NO = option_tree[i,j] = exp(-(r * h)) * (p * option_tree[i+1,j+1] + (1-p) * option_tree[i+1,j])
            EX = pmax(stock_tree[i,j]-K, 0)
            option_tree[i,j] = max(NO, EX)
         }
      }
   }
   else if(type == "put"){
      #Terminal nodes
      option_tree[nrow(option_tree),] = pmax(K-stock_tree[nrow(stock_tree),], 0)
      
      #Backward operations
      for(i in (nrow(option_tree)-1):1){
         for(j in 1:i){
            NO = option_tree[i,j] = exp(-(r * h)) * (p * option_tree[i+1,j+1] + (1-p) * option_tree[i+1,j])
            EX = pmax(K-stock_tree[i,j], 0)
            option_tree[i,j] = max(NO, EX)
         }
      }
   }
   option_tree
}

RP_stock_tree_A = function(S, K, N, h, r, delta, sigma, type){
   u = exp((r - delta)*h + sigma*sqrt(h))
   d = exp((r - delta)*h - sigma*sqrt(h))
   stock_tree = stock_tree(S, N, u, d)
   option_tree = American_option(S, K, N, h, r, delta, sigma, type)
   
   RP_stock_tree = matrix(0, nrow=nrow(stock_tree), ncol=ncol(stock_tree))
   
   for(i in (nrow(option_tree)-1):1){
      for(j in 1:i){
         RP_stock_tree[i,j] = exp(-(delta * h)) * (option_tree[i+1,j+1] - option_tree[i+1,j]) / (stock_tree[i+1,j+1] - stock_tree[i+1,j])
      }
   }
   RP_stock_tree
}

RP_bond_tree_A = function(S, K, N, h, r, delta, sigma, type){
   u = exp((r - delta)*h + sigma*sqrt(h))
   d = exp((r - delta)*h - sigma*sqrt(h))
   stock_tree = stock_tree(S, N, u, d)
   option_tree = American_option(S, K, N, h, r, delta, sigma, type)
   
   RP_bond_tree = matrix(0, nrow=nrow(stock_tree), ncol=ncol(stock_tree))
   
   for(i in (nrow(option_tree)-1):1){
      for(j in 1:i){
         RP_bond_tree[i,j] = exp(-(r * h)) * (option_tree[i+1,j] * u - option_tree[i+1,j+1] * d) / (u - d)
      }
   }
   RP_bond_tree
}

American put price would be:

American_put = American_option(10, 10, 250, 1/365, 0.01, 0, 0.15, "put")
American_put[1,1]
## [1] 0.4653426

And American call price would be:

American_call = American_option(10, 10, 250, 1/365, 0.01, 0, 0.15, "call")
American_call[1,1]
## [1] 0.5285973

Discrete Dividends

A very simple way to incorporate dividends is to assume a constant dividend yield at the payment date. A feature of this model is that the tree for the ex-dividend stock price \(S_t\) is recombining. Everything works out the same as before, except for \(u\) and \(d\), which now are defined by \(u= 1/d = e^{\sigma\sqrt{h}}\), and \(p^* = \frac{e^{rh}-d}{u-d}\)._**

Using a binomial model in which the ex-dividend price of the stock evolves according to: \[S_{t+1} = (1-I_{\{t+1\in \mathbb{D}\}}\delta)S_t \xi_{t+1}\] where \(\xi\) is a constant dividend yield, \(\mathbb{D}\subseteq\{1,...,T\}\) is the set of dividend dates, and \(\xi_{t+1}\in\{u,d\}\). The variable \(I_{\{t+1\in \mathbb{D}\}}\) takes the value of 1 if \(t+ 1\) is a dividend date and 0 otherwise. An important feature of this model is that the tree for the ex-dividend stock price is recombining.

R code which takes as inputs:

  • The payoff function \(g(S)\)

  • The nature of the derivative, i.e. European or American

  • The number of periods \(N\)

  • The interest rate \(r\)

  • The length of the period \(h\)

  • The up and down factors \(u\) and \(d\)

  • The set of dividend dates \(\mathbb{D}\subseteq\{1,...,N\}\)

  • The dividend yield \(\delta\)

  • The initial ex-dividend stock price \(S_0\)

and which outputs the price as well as the composition of the replicating portfolio at every node of the tree and which also determines the optimal exercise dates.

stock_tree_div = function(S, N, u, d, delta, D){
   tree = matrix(0, nrow=N+1, ncol=N+1)
   tree[1,1] = S
   
   for(i in 2:(N+1)){
      for(j in 1:i){
         if(j != i) {tree[i,j] = tree[i-1,j] * d}
         else if(j == i) {tree[i,j] = tree[i-1,j-1] * u}
            
         if((i-1) %in% D){
            tree[i,j] = tree[i,j] * (1-delta)
         }
      }
   }
   tree
}

Div_option = function(f, S, K, N, h, r, D, delta, sigma, type){
   u = exp(sigma*sqrt(h))
   d = 1/u
   stock_tree = stock_tree_div(S, N, u, d, delta, D)
   p = p(r, delta, h, u, d)
   option_tree = matrix(0, nrow=nrow(stock_tree), ncol=ncol(stock_tree))
   
   if(type == "European"){
      #Terminal nodes
      option_tree[nrow(option_tree),] = unlist(lapply(stock_tree[nrow(stock_tree),], f, K=K))
      
      #Backward operations
      for(i in (nrow(option_tree)-1):1){
         for(j in 1:i){
            option_tree[i,j] = exp(-(r * h)) * (p * option_tree[i+1,j+1] + (1-p) * option_tree[i+1,j])
         }
      }
   }
   
   else if(type == "American"){
      #Terminal nodes
      option_tree[nrow(option_tree),] = unlist(lapply(stock_tree[nrow(stock_tree),], f, K=K))
      
      #Backward operations
      for(i in (nrow(option_tree)-1):1){
         for(j in 1:i){
            NO = option_tree[i,j] = exp(-(r * h)) * (p * option_tree[i+1,j+1] + (1-p) * option_tree[i+1,j])
            EX = f(stock_tree[i,j], K)
            option_tree[i,j] = max(NO, EX)
         }
      }
   }
   option_tree
}

#test
#stock_tree_div(4, 3,  u = exp(0.2*sqrt(0.25)), d = 1/exp(0.2*sqrt(0.25)), 0.05, c(1,3))

#Div_option = function(f, S, K, N, h, r, D, delta, sigma, type)
#Div_option(call_payoff, 4, 4, 3, 0.25, 0.02, c(1, 3), 0.05, 0.2, "American")

#a = p(0.02, 0.05, 0.25, exp(0.2*sqrt(0.25)), 1/exp(0.2*sqrt(0.25)))
#exp(-0.02*0.25) * (a*1.87299+(1-a)*0.989667 )

RP_stock_tree_div = function(f, S, K, N, h, r, D, delta, sigma, type){
   u = exp(sigma*sqrt(h))
   d = 1/u
   stock_tree = stock_tree_div(S, N, u, d, delta, D)
   option_tree = Div_option(f, S, K, N, h, r, D, delta, sigma, type)
   
   RP_stock_tree = matrix(0, nrow=nrow(stock_tree), ncol=ncol(stock_tree))
   
   for(i in (nrow(option_tree)-1):1){
      for(j in 1:i){
         RP_stock_tree[i,j] = exp(-(delta * h)) * (option_tree[i+1,j+1] - option_tree[i+1,j]) / (stock_tree[i+1,j+1] - stock_tree[i+1,j])
      }
   }
   RP_stock_tree
}

RP_bond_tree_div = function(S, K, N, h, r, delta, sigma, type){
   u = exp(sigma*sqrt(h))
   d = 1/u
   stock_tree = stock_tree_div(S, N, u, d, delta, D)
   option_tree = Div_option(f, S, K, N, h, r, D, delta, sigma, type)
   
   RP_bond_tree = matrix(0, nrow=nrow(stock_tree), ncol=ncol(stock_tree))
   
   for(i in (nrow(option_tree)-1):1){
      for(j in 1:i){
         RP_bond_tree[i,j] = exp(-(r * h)) * (option_tree[i+1,j] * u - option_tree[i+1,j+1] * d) / (u - d)
      }
   }
   RP_bond_tree
}

Early_exercise_day = function(f, S, K, N, h, r, D, delta, sigma){
   u = exp(sigma*sqrt(h))
   d = 1/u
   stock_tree = stock_tree_div(S, N, u, d, delta, D)
   p = p(r, delta, h, u, d)
   option_tree = matrix(0, nrow=nrow(stock_tree), ncol=ncol(stock_tree))
   early_node = matrix(0, nrow=nrow(stock_tree), ncol=ncol(stock_tree))
   
   #Terminal nodes
   option_tree[nrow(option_tree),] = unlist(lapply(stock_tree[nrow(stock_tree),], f, K=K))
   
   #Backward operations
   for(i in (nrow(option_tree)-1):1){
      for(j in 1:i){
         NO = option_tree[i,j] = exp(-(r * h)) * (p * option_tree[i+1,j+1] + (1-p) * option_tree[i+1,j])
         EX = f(stock_tree[i,j], K)
         if(EX > NO){
            early_node[i,j] = 1
         }
      }
   }
   early_node
}

#test
#Early_exercise_day(call_payoff, 4, 4, 3, 0.25, 0.02, c(1, 3), 0.05, 0.2)

Applying your code to compute the initial value of an American put and an American call with strike \(K=10\) when \(r = 0.02, u = 1/d = e^{0.2\sqrt{h}}, h = 1/365, S_0 = 10, N=200\) periods, \(\delta=0.05\), and \(\mathbb{D}=\{50,100,150\}\).

#Put option payoff function
put_payoff = function(l, K){
   payoff = c()
   for (i in 1:length(l)){
      if (l[i] < K) {payoff[i] = K - l[i]}
      else {payoff[i] = 0}
   }
   payoff
}

#Div_option = function(f, S, K, N, h, r, D, delta, sigma, type)
American_put = Div_option(put_payoff, 10, 10, 200, 1/365, 0.02, c(50,100,150), 0.05, 0.2, "American")
American_put[1,1]
## [1] 1.632429
#Call option payoff function
call_payoff = function(l, K){
   payoff = c()
   for (i in 1:length(l)){
      if (l[i] > K) {payoff[i] = l[i] - K}
      else {payoff[i] = 0}
   }
   payoff
}

#Div_option = function(f, S, K, N, h, r, D, delta, sigma, type)
American_call = Div_option(call_payoff, 10, 10, 200, 1/365, 0.02, c(50,100,150), 0.05, 0.2, "American")
American_call[1,1]
## [1] 0.3042862

So the initial price of American put would be 1.632429, and the initial price of American call would be 0.3042862.

Applying code to compute the initial price of an American straddle with the same parameters as above.

#Straddle option payoff function
straddle_payoff = function(l, K){
   payoff = c()
   for (i in 1:length(l)){
      if (l[i] > K) {payoff[i] = l[i] - K}
      else {payoff[i] = K - l[i]}
   }
   payoff
}

#Div_option = function(f, S, K, N, h, r, D, delta, sigma, type)
American_straddle = Div_option(straddle_payoff, 10, 10, 200, 1/365, 0.02, c(50,100,150), 0.05, 0.2, "American")
American_straddle[1,1]
## [1] 1.772746

So the initial price of American straddle would be \(1.772746 < 1.632429 + 0.3042862\). The straddle price is smaller than the sum of call and put because American options can be exercised individually, but once the straddle is exercised, the entire package of synthetic put and call is terminated. Therefore, in an early exercise of straddle, you loss the time value of information in the options.

Pricing Asian Options by Monte Carlo

An Asian option is an option on the time average of the underlying asset option with discrete arithmetic averaging. An option with maturity of \(T\) years and strike \(K\) has the payoff \[max \Big\{\frac{1}{N} \sum_{i=1}^{N}S(t_i)-K,0\Big\}\] where \(t_i=i\times h and h=T\mathcal{N}\). Furthermore, assume \(S(t_0)=200, r=0.02, \sigma=0.20, K=220, T=1,\) and \(N=365\).

Pricing the option by Monte Carlo simulation using 100,000 paths.

Using the following stock price distribution assumption to run Monte Carlo simulation: \[\ln{S_t} \stackrel{}{\sim} N( \ln{S_0} + (r - \delta - \sigma ^2) T, \sigma^2 T)\]

stock_price = 200
strike_price = 220
r = 0.02
dividend = 0
sigma = 0.2
Time = 1
N = 365
num_sim = 100000
num_rep = 20

h = Time/N
u = exp((r-dividend)*h + sigma * (h ^ 0.5))
d = exp((r - dividend)*h - sigma * (h ^ 0.5))
prob = (exp((r - dividend) * h) - d) / (u - d)

#lnSt ~ N(lnSo + (r-div - sigma^2 / 2)T, sigma^2 T)
mu =  (r - dividend - ((sigma) ** 2 / 2)) * h
std = sigma * (h ** 0.5)

Time_0 = rep(log(stock_price),num_sim)

#Parallel computing function
sim_one_node = function(last_node){
  one_node = rnorm(num_sim,mu,std)
  new_node = last_node + one_node
  return(new_node)
}  
  
#Monte Carlo Simulation function
sim_one_matrix = function(){
  last = Time_0
  ms = as.data.frame(matrix(numeric(0),ncol=1,nrow = num_sim))
  ms[1] = Time_0
  for (i in 1:N){   #run 365 times
    out = sim_one_node(last)
    ms[i+1] = out
    last = out
  }
  ms = exp(ms)  #turn back from lnSt to St
  avg = rowMeans(ms) # take average of each path
  payoff_T = pmax((avg - strike_price), 0)  #max(Avg - K, 0)
  payoff_0 = payoff_T * exp(-r * Time)  #discount back
  avg_payoff = mean(payoff_0)  #The average price of 100000 path
  return(avg_payoff)
}
sim_one_matrix()
## [1] 3.253053

The average of 20 Monte Carlo simulations to get the 95% confidence interval. Confidence interval = mean \(\pm\) 2 std

payoffs = c()
for (i in 1:num_rep){
  payoffs[i] = sim_one_matrix()
}

total_mean = mean(payoffs)
total_sd = sd(payoffs)
interval_down = total_mean - 2 * total_sd
interval_up = total_mean + 2 * total_sd
interval_down
## [1] 3.207652
interval_up
## [1] 3.318147

Pricing American Options using the Longstaff-Schwartz LeastSquare Method

Assuming \(S_0 = 200, r = 0.1, and \sigma = 0.3\). Building a code to price an American put option with \(K = 220\) and \(T = 1\). The code takes as inputs the number of steps \(N\) and the number of simulated paths.

Pricing the option with \(N = 250\) and 100,000 paths.

#Step1 build up simulation stocks
#Monte Carlo Simulation
#lnSt ~ N(lnSo + (r-div - sigma^2 / 2)T, sigma^2 T)

set_up = function(r,dividend,sigma,Time,N){
  h = Time/N
  mu =  (r - dividend - ((sigma) ** 2 / 2)) * h
  std = sigma * (h ** 0.5)
  return(c(h,mu,std))
}

#Parallel computing function
sim_one_node = function(num_sim,mu,std,last_node){
  one_node = rnorm(num_sim,mu,std)
  new_node = last_node + one_node
  return(new_node)
}  

#simulate all paths
sim_one_matrix = function(setup,S,K,N,num_sim){
  h = setup[1]
  mu = setup[2]
  std = setup[3]
  Time_0 = rep(log(S),num_sim)
  last = Time_0
  prices = as.data.frame(matrix(numeric(0),ncol=1,nrow = num_sim))
  prices[1] = Time_0
  for (i in 1:N){   #run 365 times
    out = sim_one_node(num_sim,mu,std,last)
    prices[i+1] = out
    last = out
  }
  prices = exp(prices)  #turn back from lnSt to St
  return(prices)
}


# using Longstaff and Schwartz method
backward = function(prices_df,K,N,num_sim,r,setup){
  h = setup[1]
  one_discount = exp(-r * h)
  
  cash_flow_df = as.data.frame(matrix(numeric(0),ncol=1,nrow = num_sim))
  last_t_n = t_n = pmax(K - prices_df[,ncol(prices_df)], 0)
  cash_flow_df[1] = last_t_n
  
  #Find X,Y and do regression
  for(i in 1 : (N-1)){
    t_n = cash_flow_df[1]
    X = prices_df[,ncol(prices_df) - i]
    X[X>K] = NA  #It's a put, only cares about X<K
    EX = K- X
    Y = as.numeric(unlist(t_n * one_discount))
    X_square = X ^ 2
    out = lm(Y~ (X + X_square), na.action = na.omit) #run regression
    NO = as.numeric(out$coef[1]) + X * as.numeric(out$coef[2]) + X_square * as.numeric(out$coef[3])
    t_n_minus = ifelse(EX > NO, EX,0) #if Exercise > Not exercise, t-1 = exercise, otherwise = 0
    t_n_minus[is.na(t_n_minus)] = 0
    cash_flow_df = cbind(t_n_minus,cash_flow_df)
  }
  cash_flow_df[cash_flow_df == 0] = NA
  cash_flow_df = cash_flow_df[rowSums(is.na(cash_flow_df)) != ncol(cash_flow_df), ]# delete row with all NA
  payoff_sum = 0
  for (i in 1: nrow(cash_flow_df)){
    EX_time = min(which(!is.na(cash_flow_df[i,])))  #The earliest exercise time for each path
    payoff_sum = payoff_sum + cash_flow_df[i,EX_time] * one_discount^EX_time  #Discount the EX value back
  }
  price_0 = payoff_sum/num_sim
  return(price_0)
}

setup = set_up(0.1,0,0.3,1,250)
prices_df = sim_one_matrix(setup,200,220,250,100000)
put_0 = backward(prices_df,220,250,100000,0.1,setup)
put_0
## [1] 19.90479

Keeping \(N = 250\) steps and pricing the option when the number of paths takes the values 10, 100, 1,000, 10,000, and 100,000.

## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

Pricing the option when the number of steps \(N\) takes the values 3, 10, 100, 250, and 1,000.

plot(c(f,g,h,i,j),xaxt = "n",ylab = "Price", xlab = "N", main = "Amrican Put Price under different path simulation", pch = 16, col = "blue")
axis(side = 1,at = 1:5, labels = c(3,10,100,250,1000))