Nested Fors

Examples by Prof. Jacob Escobar, EGADE Business School, Tec de Monterrey

Typical convention: [i, j] -> i for rows, j for columns

A1 <- matrix(0, 4, 5)
iter <- 0
for (i in 1:4) {
  #print(paste("Row i:", i))
  for (j in 1:5)   {
    #print(paste("Colum j:", j))
    iter <- iter + 1
    A1[i, j] <- iter
    #print(A1)
  }
}
A1

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    2    3    4    5
## [2,]    6    7    8    9   10
## [3,]   11   12   13   14   15
## [4,]   16   17   18   19   20

Swapping innner and outer loops:

A2 <- matrix(0, 4, 5)
iter <- 0
for (j in 1:5) {
  for (i in 1:4) {
    iter <- iter + 1
    A2[i, j] <- iter
  }
}
A2

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    5    9   13   17
## [2,]    2    6   10   14   18
## [3,]    3    7   11   15   19
## [4,]    4    8   12   16   20

Referencing the diagonal, upper o lower half of the matrix * How can you write "d"s in the diagonal? * How can you write "u"s in the upper half and "l"s in the lower half?

B1 <- matrix(0, 4, 4)
for (j in 1:4) {
  for (i in 1:4) {
    if(i == j) B1[i,j] <- "d"
    if(j > i) B1[i,j] <- "u"
    if(j < i) B1[i,j] <- "l"
  }
}
B1

##      [,1] [,2] [,3] [,4]
## [1,] "d"  "u"  "u"  "u" 
## [2,] "l"  "d"  "u"  "u" 
## [3,] "l"  "l"  "d"  "u" 
## [4,] "l"  "l"  "l"  "d"

Filling the upper half only:

B2 <- matrix(NA, 4, 4)
for (j in 1:4) {
  for (i in 1:4) {
    if(j >= i) B2[i,j] <- 0
  }
}
B2

##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    0
## [2,]   NA    0    0    0
## [3,]   NA   NA    0    0
## [4,]   NA   NA   NA    0

Alternatives to filling in only the upper half of the matrix (without the ifs):

C1 <- matrix(NA, 4, 4)
for (j in 1:4) {
  for (i in 1:j) { # Notice 1:j
    C1[i,j] <- 0
  }
}
C1

##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    0
## [2,]   NA    0    0    0
## [3,]   NA   NA    0    0
## [4,]   NA   NA   NA    0

C2 <- matrix(NA, 4, 4)
for (i in 1:4) {
  for (j in i:4) { # Notice i:4
    C2[i,j] <- 0
  }
}
C2

##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    0
## [2,]   NA    0    0    0
## [3,]   NA   NA    0    0
## [4,]   NA   NA   NA    0

Using vectors to fill rows entirely:

D1 <- matrix(NA, 4, 4)
for (i in 1:4) {
  D1[i, ] <- letters[1:4] # fill entire row
}
D1

##      [,1] [,2] [,3] [,4]
## [1,] "a"  "b"  "c"  "d" 
## [2,] "a"  "b"  "c"  "d" 
## [3,] "a"  "b"  "c"  "d" 
## [4,] "a"  "b"  "c"  "d"

Using vectors to fill columns entirely:

D2 <- matrix(NA, 4, 4)
for (j in 1:4) {
  D2[, j] <- letters[1:4] # fill entire column
}
D2

##      [,1] [,2] [,3] [,4]
## [1,] "a"  "a"  "a"  "a" 
## [2,] "b"  "b"  "b"  "b" 
## [3,] "c"  "c"  "c"  "c" 
## [4,] "d"  "d"  "d"  "d"

Binomial Trees

Spot Prices Simulation

Example with 3 steps (\(N = 3\)):

			\(S_0u^3\)
		\(S_0u^2\)
	\(S_0u\)		\(S_0u^2d\)
\(S_0\)		\(S_0ud\)
	\(S_0d\)		\(S_0ud^2\)
		\(S_0d^2\)
			\(S_0d^3\)

where:

\(S_0\) = spot price
\(u\) = factor by which \(S_0\) goes up in \(\Delta t\)
\(d\) = factor by which \(S_0\) goes down in \(\Delta t\)
\(T\) = time to maturity
\(N\) = number of steps

\[ \Delta t = \frac{T}{N} \]

\[ u = e^{\sigma\sqrt{\Delta t}} \]

\[ d = e^{-\sigma\sqrt{\Delta t}} = \frac{1}{u}\]

Validation tip: 1. What range of values can \(u\) have? 2. What range of values can \(d\) have?

Inputs: * n : number of steps for valuation (quick tests n = 10, suggested n = 1000) * s : underlying asset’s spot price * t : option’s time to maturity in years (e.g. 3 months = 3/12 years) * v : underlying asset’s returns annualized volatility

t <- 0.50
n <- 10
dt <- _____
dt

## [1] 0.05

s <- 52
v <- 0.20
u <- _____
u

## [1] 1.15191

d <- _____
d

## [1] 0.8681234

Is \(u > 1\)?
Is \(0 < d < 1\)?

Generation of underlier’s prices tree:

st matrix initialization, how many rows and columns?

st <- matrix(NA, _____, _____)
st

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [2,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [3,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [4,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [5,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [6,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [7,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [8,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [9,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
## [10,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
## [11,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA

Remember that we only want the top half-matrix.

for(j in 1:(_____)) {
  for (i in 1:_____) {
    st[i, j] <- s*u^(_____)*d^(_____)
  }
}
st

##       [,1]     [,2]     [,3]     [,4]     [,5]      [,6]      [,7]      [,8]
##  [1,]   52 59.89932 68.99861 79.48019 91.55402 105.46198 121.48270 139.93712
##  [2,]   NA 45.14242 52.00000 59.89932 68.99861  79.48019  91.55402 105.46198
##  [3,]   NA       NA 39.18919 45.14242 52.00000  59.89932  68.99861  79.48019
##  [4,]   NA       NA       NA 34.02106 39.18919  45.14242  52.00000  59.89932
##  [5,]   NA       NA       NA       NA 29.53448  34.02106  39.18919  45.14242
##  [6,]   NA       NA       NA       NA       NA  25.63957  29.53448  34.02106
##  [7,]   NA       NA       NA       NA       NA        NA  22.25831  25.63957
##  [8,]   NA       NA       NA       NA       NA        NA        NA  19.32296
##  [9,]   NA       NA       NA       NA       NA        NA        NA        NA
## [10,]   NA       NA       NA       NA       NA        NA        NA        NA
## [11,]   NA       NA       NA       NA       NA        NA        NA        NA
##            [,9]     [,10]     [,11]
##  [1,] 161.19496 185.68207 213.88902
##  [2,] 121.48270 139.93712 161.19496
##  [3,]  91.55402 105.46198 121.48270
##  [4,]  68.99861  79.48019  91.55402
##  [5,]  52.00000  59.89932  68.99861
##  [6,]  39.18919  45.14242  52.00000
##  [7,]  29.53448  34.02106  39.18919
##  [8,]  22.25831  25.63957  29.53448
##  [9,]  16.77472  19.32296  22.25831
## [10,]        NA  14.56253  16.77472
## [11,]        NA        NA  12.64207

Validation tips: 1. Do the prices go up and down as expected? 2. Do you see the recombinant feature?

Option Prices Backwards Induction

If he risk-neutral probability of u is p and the risk-neutral probability of d is q, then:

\[ p = \frac{e^{r\Delta t}-d}{u - d}\]

\[ q = \frac{u - e^{r\Delta t}}{u - d}\]

Validation tips:

What range of values can p and q have?
What condition do they have to meet?

Additional Inputs: * type : call (c) or put (p); * style : european (e) or american (a) * k : option’s strike price * r : annualized risk-free rate

k <- 50
r <- 0.10
p <- (_____ - d)/(_____ - d)
p

## [1] 0.6453713

q <- _____
q

## [1] 0.3546287

Is \(0 < p < 1\)?
Is \(0 < q < 1\)?
Is \(p + q = 1\)?

Generation of option’s prices tree:

type <- "c"
style <- "e"

op matrix initialization, how many rows and columns?

op <- matrix(NA, _____, _____)
op

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [2,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [3,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [4,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [5,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [6,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [7,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [8,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
##  [9,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
## [10,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA
## [11,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA

Last column: calculation of each scenario’s payout at T.

for(i in 1:(_____)) {
  if(type == "c") op[i, _____] <- max(st[i, _____] - k, 0)
  if(type == "p") op[i, _____] <- max(k - st[i, _____], 0)
}
op

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]     [,11]
##  [1,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA 163.88902
##  [2,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA 111.19496
##  [3,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA  71.48270
##  [4,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA  41.55402
##  [5,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA  18.99861
##  [6,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA   2.00000
##  [7,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA   0.00000
##  [8,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA   0.00000
##  [9,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA   0.00000
## [10,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA   0.00000
## [11,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA   0.00000

Rest of columns: backwards induction discounting the payout to the present-time (European option). If American option, also test for early-exercise.

for(j in _____:1) {
  for (i in 1:j) {
    op[i, j] <- (op[_____]*p + op[_____]*q)*exp(-r*dt)
    if(style == "a" & type == "c") op[i, j] <- max(_____, st[_____] - k)    
    if(style == "a" & type == "p") op[i, j] <- max(_____, k - st[_____])    
  }
}
op

##           [,1]     [,2]     [,3]      [,4]      [,5]       [,6]       [,7]
##  [1,] 35.32178 42.71976 51.17749 60.735938 71.431098 83.3069365 96.4312627
##  [2,]       NA 22.35781 27.93177 34.505903 42.130815 50.8284882 60.6001247
##  [3,]       NA       NA 12.53006 16.362626 21.117425 26.8978462 33.7640037
##  [4,]       NA       NA       NA  5.732462  7.940879 10.8964047 14.7826517
##  [5,]       NA       NA       NA        NA  1.794498  2.6745031  3.9780299
##  [6,]       NA       NA       NA        NA        NA  0.2183841  0.3400813
##  [7,]       NA       NA       NA        NA        NA         NA  0.0000000
##  [8,]       NA       NA       NA        NA        NA         NA         NA
##  [9,]       NA       NA       NA        NA        NA         NA         NA
## [10,]       NA       NA       NA        NA        NA         NA         NA
## [11,]       NA       NA       NA        NA        NA         NA         NA
##              [,8]        [,9]      [,10]     [,11]
##  [1,] 110.9075897 126.8728890 144.477997 163.88902
##  [2,]  71.4495174  83.4207533  96.627507 111.19496
##  [3,]  41.7124022  50.6735589  60.565544  71.48270
##  [4,]  19.7763417  25.9940184  33.387867  41.55402
##  [5,]   5.9038259   8.7406262  12.905729  18.99861
##  [6,]   0.5295959   0.8247198   1.284305   2.00000
##  [7,]   0.0000000   0.0000000   0.000000   0.00000
##  [8,]   0.0000000   0.0000000   0.000000   0.00000
##  [9,]          NA   0.0000000   0.000000   0.00000
## [10,]          NA          NA   0.000000   0.00000
## [11,]          NA          NA         NA   0.00000

Validation tips: 1. Does the value on the top-left corner make sense? 2. Do you see 0’s for part of the tree, and do they make sense based on the type of option (call vs put)?

Custom Function

bino <- function(n = 2000, type = "c", style = "e",
                 s = 52, k = 50, t = 0.5, v = 0.20, r = 0.10) {
  
  dt <- t/n
  u <- exp(v*sqrt(dt))
  d <- 1/u
  p <- (exp(r*dt) - d)/(u - d)
  q <- 1 - p

  # generation of underlier's prices tree
  st <- matrix(NA, n+1, n+1)
  for(j in 1:(n+1)) {
    for (i in 1:j) {
      st[i, j] <- s*u^(j-i)*d^(i-1)
    }
  }
  
  # generation of option prices tree
  # last column: calculation of each scenario's payout at T
  op <- matrix(NA, n+1, n+1)
  # notice how we now fill the entire column instead of using a for loop
  if (type == "c") op[, n+1] <- pmax(st[, n+1] - k, 0)
  if (type == "p") op[, n+1] <- pmax(k - st[, n+1], 0)

  # rest of columns: backwards induction discounting the payout
  # to the present-time, while testing for early-exercise
  for(j in n:1) {
    for (i in 1:j) {
      op[i, j] <- (op[i, j+1]*p + op[i+1, j+1]*q)*exp(-r*dt)
      if(style == "a" & type == "c") op[i, j] <- max(op[i, j], st[i, j] - k)    
      if(style == "a" & type == "p") op[i, j] <- max(op[i, j], k - st[i, j])    
    }
  }
  
  # How do we return the final option price only,
  # and not the whole option prices tree?
  return(___)
}

Verify that the function in the workspace:

bino

## function(n = 2000, type = "c", style = "e",
##                  s = 52, k = 50, t = 0.5, v = 0.20, r = 0.10) {
##   
##   dt <- t/n
##   u <- exp(v*sqrt(dt))
##   d <- 1/u
##   p <- (exp(r*dt) - d)/(u - d)
##   q <- 1 - p
##   print(paste("u:", u, "d:", d, "p:", p, "q:", q))
## 
##   # generation of underlier's prices tree
##   st <- matrix(NA, n+1, n+1)
##   for(j in 1:(n+1)) {
##     for (i in 1:j) {
##       st[i, j] <- s*u^(j-i)*d^(i-1)
##     }
##   }
##   
##   # generation of option prices tree
##   # last column: calculation of each scenario's payout at T
##   op <- matrix(NA, n+1, n+1)
##   # notice how we now fill the entire column instead of using a for loop
##   if (type == "c") op[, n+1] <- pmax(st[, n+1] - k, 0)
##   if (type == "p") op[, n+1] <- pmax(k - st[, n+1], 0)
## 
##   # rest of columns: backwards induction discounting the payout
##   # to the present-time, while testing for early-exercise
##   for(j in n:1) {
##     for (i in 1:j) {
##       op[i, j] <- (op[i, j+1]*p + op[i+1, j+1]*q)*exp(-r*dt)
##       if(style == "a" & type == "c") op[i, j] <- max(op[i, j], st[i, j] - k)    
##       if(style == "a" & type == "p") op[i, j] <- max(op[i, j], k - st[i, j])    
##     }
##   }
##   
##   # How do we return the final option price only,
##   # and not the whole option prices tree?
##   return(op[1,1])
## }

Test the function with the default values:

bino()

## [1] "u: 1.0031672829348 d: 0.996842717073533 p: 0.503162321141831 q: 0.496837678858169"

## [1] 5.564866

Test the function with other values:

bino(type = "p")

## [1] "u: 1.0031672829348 d: 0.996842717073533 p: 0.503162321141831 q: 0.496837678858169"

## [1] 1.126337

Test the function with other values:

bino(type = "p", style = "a")

## [1] "u: 1.0031672829348 d: 0.996842717073533 p: 0.503162321141831 q: 0.496837678858169"

## [1] 1.273054

04b Nested Fors & Binomial Trees

Jacob Escobar

31/08/2021

Nested Fors

Binomial Trees

Spot Prices Simulation

Option Prices Backwards Induction

Custom Function