Attribute VB_Name = "Ch6"
Function Black_Scholes_Call(ByVal S0 As Double, ByVal T As Double, ByVal K As Double, ByVal r As Double, ByVal delta As Double, ByVal sigma As Double) As Double
  'S0 is the initial stock price
  'T is the time to maturity
  'K is the strike price
  'r is the continuous interest rate
  'delta is the continuous dividend rate
  'sigma is the stock volatility
  Black_Scholes_Call = S0 * Exp(-delta * T) * WorksheetFunction.NormSDist(Black_Scholes_D1(S0, T, K, r, delta, sigma)) - K * Exp(-r * T) * WorksheetFunction.NormSDist(Black_Scholes_D2(S0, T, K, r, delta, sigma))
End Function

Function Black_Scholes_Put(ByVal S0 As Double, ByVal T As Double, ByVal K As Double, ByVal r As Double, ByVal delta As Double, ByVal sigma As Double) As Double
  'S0 is the initial stock price
  'T is the time to maturity
  'K is the strike price
  'r is the continuous interest rate
  'delta is the continuous dividend rate
  'sigma is the stock volatility
  Black_Scholes_Put = K * Exp(-r * T) * WorksheetFunction.NormSDist(-Black_Scholes_D2(S0, T, K, r, delta, sigma)) - S0 * Exp(-delta * T) * WorksheetFunction.NormSDist(-Black_Scholes_D1(S0, T, K, r, delta, sigma))
End Function

Private Function Black_Scholes_D1(ByVal S0 As Double, ByVal T As Double, ByVal K As Double, ByVal r As Double, ByVal delta As Double, ByVal sigma As Double) As Double
  Black_Scholes_D1 = (Log(S0 / K) + (r - delta + sigma ^ 2 / 2) * T) / sigma / Sqr(T)
End Function

Private Function Black_Scholes_D2(ByVal S0 As Double, ByVal T As Double, ByVal K As Double, ByVal r As Double, ByVal delta As Double, ByVal sigma As Double) As Double
  Black_Scholes_D2 = Black_Scholes_D1(S0, T, K, r, delta, sigma) - sigma * Sqr(T)
End Function

Function Black_Scholes_Call_Implied_Vol(ByVal S0 As Double, ByVal T As Double, ByVal K As Double, ByVal r As Double, ByVal delta As Double, ByVal price As Double) As Double
  Dim param(5)
  param(0) = S0
  param(1) = T
  param(2) = K
  param(3) = r
  param(4) = delta
  Black_Scholes_Call_Implied_Vol = invFuncS(price, "Black_Scholes_Call_2", , 0.01, , param)
End Function

Private Function Black_Scholes_Call_2(ByVal sigma As Double, param()) As Double
  Black_Scholes_Call_2 = Black_Scholes_Call(param(0), param(1), param(2), param(3), param(4), sigma)
End Function

Function invFunc(ByVal y As Double, ByVal f As String, Optional ByVal sigFig As Integer = 6, _
    Optional ByVal initGuessValue As Double = 1, Optional param) As Double
  'Newton's method to approximate x from 0 = f(x) - y
  If sigFig <= 0 Then
    sigFig = 6
  ElseIf sigFig > 15 Then
    sigFig = 15
  End If
  
  Dim oldX As Double, newX As Double, df As Double
  oldX = initGuessValue
  Do
    If IsArray(param) Then
      df = Application.Run("derivative", f, oldX, param)
    Else
      df = Application.Run("derivative", f, oldX)
    End If
    
    If df = 0 Then Exit Do
    
    If IsArray(param) Then
      newX = oldX - (Application.Run(f, oldX, param) - y) / df
    Else
      newX = oldX - (Application.Run(f, oldX) - y) / df
    End If
    
    If magnitude(newX) < -100 Then
      invFunc = 0
      Exit Function
    End If
    
    Dim compX As Double
    compX = Round(newX * 10# ^ (sigFig - magnitude(newX)))
    
    If compX = Round(oldX * 10# ^ (sigFig - magnitude(newX))) Then
      invFunc = compX / 10# ^ (sigFig - magnitude(newX))
      Exit Function
    End If
    
    oldX = newX
  Loop
  invFunc = Round(newX * 10# ^ (sigFig - magnitude(newX))) / _
      10# ^ (sigFig - magnitude(newX))
End Function

Function magnitude(ByVal x As Double) As Integer
  'the order of magnitude of a number
  x = Abs(x)
  If x = 0 Then
    magnitude = 0
  Else
    magnitude = Int(Log(x) / Log(10#))
  End If
End Function

Function derivative(ByVal f As String, ByVal x As Double, Optional param) As Double
  'Approximate the derivative of f at point x
  Dim dx As Double
  If x = 0 Then dx = 2# ^ -30 Else dx = x / (2# ^ 30)
  
  If IsArray(param) Then
    derivative = (Application.Run(f, x + dx, param) - Application.Run(f, x, param)) / dx
  Else
    derivative = (Application.Run(f, x + dx) - Application.Run(f, x)) / dx
  End If
End Function

Function invFuncS(ByVal y As Double, ByVal f As String, Optional ByVal sigFig As Integer = 6, _
    Optional ByVal initGuessValue1 As Double = 0.01, Optional ByVal initGuessValue2 As Double _
    = 1, Optional param) As Double
  'Secant method to approximate x from 0 = f(x) - y
  If sigFig <= 0 Then
    sigFig = 6
  ElseIf sigFig > 15 Then
    sigFig = 15
  End If
  
  Dim oldX As Double, olderX As Double, newX As Double
  olderX = initGuessValue1
  oldX = initGuessValue2
  
  Do
    If IsArray(param) Then
      newX = oldX - (oldX - olderX) / (Application.Run(f, oldX, param) - Application.Run _
          (f, olderX, param)) * (Application.Run(f, oldX, param) - y)
    Else
      newX = oldX - (oldX - olderX) / (Application.Run(f, oldX) - Application.Run(f, olderX)) _
          * (Application.Run(f, oldX) - y)
    End If
    
    If magnitude(newX) < -100 Then
      invFuncS = 0
      Exit Function
    End If
    
    Dim compX As Double
    compX = Round(newX * 10# ^ (sigFig - magnitude(newX)))
    
    If compX = Round(oldX * 10# ^ (sigFig - magnitude(newX))) Then
      invFuncS = compX / 10# ^ (sigFig - magnitude(newX))
      Exit Function
    End If
    olderX = oldX
    oldX = newX
  Loop
End Function

Private Function ch6_1_func(x)
  func = x ^ 2 + 2 * x + 1
End Function

Sub ch6_1()
  MsgBox invFunc(0, "ch6_1_func", 12) & vbCrLf & invFuncS(0, "ch6_1_func", 12)
End Sub

Sub ch6_2()
  Dim u As Double, strArray(9) As String, i As Integer
  Randomize

  'simulate 10 standard normal random variable
  For i = 0 To 9
    'generate a U(0,1) random variable
    u = Rnd
    'calculate the simulated value of a normal random variable
    strArray(i) = Round(WorksheetFunction.NormSInv(u), 6)
  Next i
  
  MsgBox Join(strArray, vbCrLf)
End Sub

Sub ch6_2_ex1()
  Randomize
  Dim dblArray(1 To 10000, 1 To 2) As Double
  Cells(1, 1) = "# trial"
  Cells(1, 2) = "Exp(1)"
  For i = 1 To 10000
    dblArray(i, 1) = i
    dblArray(i, 2) = WorksheetFunction.GammaInv(Rnd, 1, 1)
    DoEvents
  Next i
  Range(Cells(2, 1), Cells(10001, 2)) = dblArray
End Sub

Sub ch6_2_ex1_another_method()
  Randomize
  Cells(1, 1) = "# trial"
  Cells(1, 2) = "Exp(1)"
  Cells(2, 1) = "generating..."
  Application.Calculation = xlCalculationManual
  Application.ScreenUpdating = False
  For i = 1 To 10000
    Cells(i + 1, 1) = i
    Cells(i + 1, 2) = WorksheetFunction.GammaInv(Rnd, 1, 1)
    DoEvents
  Next i
  Application.Calculation = xlCalculationAutomatic
  Application.ScreenUpdating = True
End Sub

Sub ch6_2_ex2()
  Randomize
  Dim dblArray(1 To 5000, 1 To 2) As Double
  Cells(1, 1) = "# trial"
  Cells(1, 2) = "Exp(1)"
  For i = 1 To 5000
    dblArray(i, 1) = i
    Dim rndNum As Double
    rndNum = Rnd
    
    'notice a jump between 0.5 and 0.75 (at X = 1)
    If 0.5 < rndNum And rndNum < 0.75 Then
      dblArray(i, 2) = 1
    ElseIf rndNum <= 0.5 Then
      dblArray(i, 2) = invFuncS(rndNum, "ch6_2_ex2_func")
    Else
      dblArray(i, 2) = invFuncS(rndNum, "ch6_2_ex2_func")
    End If
    DoEvents
  Next i
  Range(Cells(2, 1), Cells(5001, 2)) = dblArray
End Sub

Private Function ch6_2_ex2_func(x)
  If x < 0 Then
    ch6_2_ex2_func = 0
  ElseIf 0 <= x And x < 1 Then
    ch6_2_ex2_func = 0.5 * x
  ElseIf 1 <= x And x < 2 Then
    ch6_2_ex2_func = 0.5 + 0.25 * x
  Else
    ch6_2_ex2_func = 1
  End If
End Function

Function discDistInv(ByVal u As Double, ByVal pmf As String, Optional param) As Long
  'determine the inverse of a discrete distribution
  Dim x As Long, oldFx As Double, newFx As Double
  x = 0
  oldFx = 0
  Do
    If IsArray(param) Then
      newFx = oldFx + Application.Run(pmf, x, param)
    Else
      newFx = oldFx + Application.Run(pmf, x)
    End If
    
    If oldFx <= u And u < newFx Then
      discDistInv = x
      Exit Function
    End If
    
    x = x + 1
    oldFx = newFx
  Loop
End Function

Private Function ch6_3_ex_func(x)
  ch6_3_ex_func = WorksheetFunction.BinomDist(x, 4, 0.6, False)
End Function

Sub ch6_3_ex()
  Randomize
  Dim intArray(1 To 10000, 1 To 2) As Integer
  Cells(1, 1) = "# trial"
  Cells(1, 2) = "B(4, 0.6)"
  For i = 1 To 10000
    intArray(i, 1) = i
    intArray(i, 2) = discDistInv(Rnd, "ch6_3_ex_func")
    DoEvents
  Next i
  Range(Cells(2, 1), Cells(10001, 2)) = intArray
End Sub

Sub ch6_3_ex_another_method()
  Randomize
  Cells(1, 1) = "# trial"
  Cells(1, 2) = "B(4, 0.6)"
  Cells(2, 1) = "generating..."
  Application.Calculation = xlCalculationManual
  Application.ScreenUpdating = False
  For i = 1 To 10000
    Cells(i + 1, 1) = i
    Cells(i + 1, 2) = discDistInv(Rnd, "ch6_3_ex_func")
    DoEvents
  Next i
  Application.Calculation = xlCalculationAutomatic
  Application.ScreenUpdating = True
End Sub

Sub ch6_4_ex()
  Randomize
  
  Dim i As Integer, j As Integer, result(1 To 10000, 1 To 2) As Double

  Cells(1, 1) = "# trial"
  Cells(1, 2) = "S"
  
  For j = 1 To 10000
    'First simulate N
    Dim n As Integer
    n = discDistInv(Rnd, "ch6_4_ex_funcN")
  
    'Generate Xi
    Dim sum As Double
    sum = 0
    For i = 1 To n
      sum = sum + invFuncS(Rnd, "ch6_4_ex_funcX")
    Next i
    
    result(j, 1) = j
    result(j, 2) = sum
    
    DoEvents
  Next j
  
  Range(Cells(2, 1), Cells(10001, 2)) = result
End Sub

Private Function ch6_4_ex_funcN(x)
  ch6_4_ex_funcN = WorksheetFunction.Poisson(x, 8, False)
End Function

Private Function ch6_4_ex_funcX(x)
  ch6_4_ex_funcX = 1 - (5000 / (x + 5000)) ^ 2
End Function

Sub ch6_5_ex1()
  Dim simValues(1 To 50000) As Double, sum As Double, i As Long
  
  'Simulate values and calculate sum
  For i = 1 To 50000
    simValues(i) = invFuncS(Rnd, "ch6_5_ex_func")
    sum = sum + simValues(i)
    DoEvents
  Next i
  
  'Calculate mean
  Dim mean As Double
  mean = sum / 50000
  
  'Calculate variance
  Dim variance As Double
  For i = 1 To 50000
    variance = variance + (simValues(i) - mean) ^ 2
    DoEvents
  Next i
  variance = variance / 49999
  
  Dim trueMean As Double, trueVar As Double
  trueMean = 250
  trueVar = 500 ^ 2 - trueMean ^ 2
  
  MsgBox "Sim Mean: " & mean & "; True Mean: " & trueMean & vbCrLf & _
         "Sim Var: " & variance & "; True Var: " & trueVar
End Sub

Private Function ch6_5_ex_func(x)
  ch6_5_ex_func = 1 - (500 / (x + 500)) ^ 3
End Function

Sub ch6_5_ex2()
  Dim i As Integer, j As Integer, total As Double
  
  '4 months
  For i = 1 To 4
    Dim sum As Double
    sum = 0
    
    For j = 1 To 10000
      sum = sum + WorksheetFunction.Max(WorksheetFunction.NormInv(Rnd, 15000, 2000) - 10000, 0)
      DoEvents
    Next j
    
    total = total + sum / 10000
  Next i
  
  MsgBox "Yearly simulated cost: " & total
End Sub

Function simBinomialStock(ByVal S0 As Double, ByVal u As Double, ByVal d As Double, _
    ByVal r As Double, Optional ByVal delta As Double = 0, Optional ByVal step As Integer = 1, _
    Optional ByVal h As Double = 1) As Double
  'S0 is the initial stock price
  'u is the ratio of the new stock price to the initial stock price if going up
  'd is the ratio of the new stock price to the initial stock price if going down
  'r is the continuous interest rate
  'delta is the continuous dividend rate
  'step is the number of steps of the binomial tree
  'h is the time period of each step
  Dim param(1), nSuccess As Integer
  
  'Get the simulated binomial variable
  param(0) = step                                 'n = step
  param(1) = (Exp((r - delta) * h) - d) / (u - d) 'p = p*
  nSuccess = discDistInv(Rnd, "simBinomialStock_func", param)
  
  'Calculate the stock price
  simBinomialStock = S0 * u ^ nSuccess * d ^ (step - nSuccess)
  
End Function

Private Function simBinomialStock_func(ByVal x As Integer, param As Variant) As Double
  'return a binomial pmf
  simBinomialStock_func = WorksheetFunction.BinomDist(x, param(0), param(1), False)
End Function

Function simBinomialStockGBM(ByVal S0 As Double, ByVal sigma As Double, _
    ByVal r As Double, Optional ByVal delta As Double = 0, Optional ByVal step As Integer = 1, _
    Optional ByVal h As Double = 1) As Double
  'S0 is the initial stock price
  'u is the ratio of the new stock price to the initial stock price if going up
  'd is the ratio of the new stock price to the initial stock price if going down
  'r is the continuous interest rate
  'delta is the continuous dividend rate
  'step is the number of steps of the binomial tree
  'h is the time period of each step
  Dim u As Double, d As Double
  u = Exp((r - delta) * h + sigma * Sqr(h))
  d = Exp((r - delta) * h - sigma * Sqr(h))
  simBinomialStockGBM = simBinomialStock(S0, u, d, r, delta, step, h)
End Function

Sub ch6_6_1_ex1()
  Randomize
  
  Dim i As Integer, j As Integer, result(1 To 10) As String
  
  For i = 1 To 10
    Dim sum As Double
    sum = 0
    For j = 1 To 5000
      sum = sum + simBinomialStock(100, 1.002, 0.998, 0.05, 0.02, 100, 0.01)
      DoEvents
    Next j
    result(i) = sum / 5000
  Next i
  
  MsgBox Join(result, vbCrLf)
End Sub

Sub ch6_6_1_ex2()
  Randomize
  
  Dim i As Integer, j As Integer, result(1 To 10) As String
  
  Dim u As Double, d As Double
  u = Exp((0.05 - 0.02) * 0.01 + 0.2 * Sqr(0.01))
  d = Exp((0.05 - 0.02) * 0.01 - 0.2 * Sqr(0.01))
  
  For i = 1 To 10
    Dim sum As Double
    sum = 0
    For j = 1 To 5000
      sum = sum + simBinomialStock(100, u, d, 0.05, 0.02, 100, 0.01)
      DoEvents
    Next j
    result(i) = sum / 5000
  Next i
  
  MsgBox Join(result, vbCrLf)
End Sub

Function simGBMStock(ByVal S0 As Double, ByVal sigma As Double, _
    ByVal r As Double, Optional ByVal delta As Double = 0, _
    Optional ByVal T As Double = 1) As Double
  'S0 is the initial stock price
  'sigma is the volatility of the stock price
  'r is the continuous interest rate
  'delta is the continuous dividend rate
  't is the elapsed time
  simGBMStock = S0 * Exp((r - delta - 0.5 * sigma ^ 2) * T + _
      sigma * WorksheetFunction.NormInv(Rnd, 0, Sqr(T)))
End Function

Sub ch6_6_2_ex()
  Randomize
  
  Dim i As Integer, j As Integer, result(1 To 10) As String
  
  
  For i = 1 To 10
    Dim sum As Double
    sum = 0
    For j = 1 To 5000
      sum = sum + simGBMStock(100, 0.2, 0.05, 0.02, 1)
      DoEvents
    Next j
    result(i) = sum / 5000
  Next i
  
  MsgBox Join(result, vbCrLf)
End Sub

Sub ch6_7_ex()
  Randomize

  Dim j As Integer
  
  Dim u As Double, d As Double
  u = Exp(0.08 * 0.01 + 0.3 * Sqr(0.01))
  d = Exp(0.08 * 0.01 - 0.3 * Sqr(0.01))
  
  Dim sum As Double, price1 As Double
  For j = 1 To 10000
    sum = sum + WorksheetFunction.Max(simBinomialStock(41, u, d, 0.08, 0, 100, 0.01) - 40, 0)
    DoEvents
  Next j
  
  price1 = sum / 10000 * Exp(-0.08)
  
  Dim price2 As Double
  sum = 0
  For j = 1 To 10000
    sum = sum + WorksheetFunction.Max(simGBMStock(41, 0.3, 0.08, 0, 1) - 40, 0)
    DoEvents
  Next j
  
  price2 = sum / 10000 * Exp(-0.08)
  
  Dim truePrice As Double
  truePrice = Black_Scholes_Call(41, 1, 40, 0.08, 0, 0.3)
  
  MsgBox "Price using Binomial Model: " & price1 & vbCrLf & _
         "Price using GBM Model: " & price2 & vbCrLf & _
         "True Price: " & truePrice
End Sub

Sub ch6_8_ex()
  Randomize
  Dim j As Integer, sum As Double
  Dim S(2) As Double
  
  For j = 1 To 10000
    S(0) = simGBMStock(41, 0.3, 0.08, 0, 1 / 3) 'First Path
    S(1) = simGBMStock(S(0), 0.3, 0.08, 0, 1 / 3) 'Second Path
    S(2) = simGBMStock(S(1), 0.3, 0.08, 0, 1 / 3)  'Third Path
    sum = sum + WorksheetFunction.Max((S(0) + S(1) + S(2)) / 3 - 40, 0)
    DoEvents
  Next j
  
  MsgBox "The simulated price is " & sum / 10000 * Exp(-0.08)
End Sub

Function Chol(cov)
  ' Cholesky decomposition
  ' Cov = L x L matrix of covariances
  Dim SumSq As Double, SumPr As Double, h As Integer
  Dim i As Integer, j As Integer, a() As Double, L As Integer
  L = UBound(cov) - LBound(cov) + 1
  ReDim a(L - 1, L - 1)

  For i = 0 To L - 1
    SumSq = 0
    For h = 0 To i - 1
      SumSq = SumSq + a(i, h) * a(i, h)
    Next h
    a(i, i) = Sqr(cov(i, i) - SumSq)
    For j = i To L - 1
      SumPr = 0
      For h = 0 To i - 1
        SumPr = SumPr + a(i, h) * a(j, h)
      Next h
      a(j, i) = (cov(i, j) - SumPr) / a(i, i)
    Next j
  Next i
  Chol = a
End Function

Function simCorrGBMStock(S0() As Double, delta() As Double, sigma() As Double, _
    corr() As Double, ByVal r As Double, Optional ByVal T As Double = 1) As Double()
  'S0 is a vector of initial stock prices
  'delta is a vector of dividend rates
  'sigma is a vector of volatilities
  'corr is the correlation matrix
  'r is the risk-free rate
  'T is the time elapsed
  Dim ST() As Double, epsillon() As Double, a() As Double, Z() As Double, n As Integer
  
  n = UBound(S0) - LBound(S0) + 1
  ReDim ST(n - 1)
  ReDim epsillon(n - 1)
  ReDim a(n - 1, n - 1)
  ReDim Z(n - 1)
  
  'Simulate n uncorrelated standard normal
  Dim K As Integer, i As Integer
  For K = 0 To n - 1
    epsillon(K) = WorksheetFunction.NormSInv(Rnd)
  Next K
  
  'Calculate the Cholesky decomposition
  a = Chol(corr)
  
  'Calculate the n correlated N(0,t)
  For K = 0 To n - 1
    For i = 0 To K
      Z(K) = a(K, i) * epsillon(i) * Sqr(T)
    Next i
  Next K
  
  'Calculate the stock prices
  For K = 0 To n - 1
    ST(K) = S0(K) * Exp((r - delta(K) - 0.5 * sigma(K) ^ 2) * T + _
        delta(K) * Z(K))
  Next K
  
  simCorrGBMStock = ST
End Function

Function European_Basket_Call_Option_MC(S0() As Double, delta() As Double, sigma() As Double, corr() As Double, w() As Double, ByVal r As Double, ByVal K As Double, ByVal T As Double) As Double
  'S0 is a vector of initial stock prices
  'delta is a vector of dividend rates
  'sigma is a vector of volatilities
  'corr is the correlation matrix
  'w is a vector of weights
  'r is the risk-free rate
  'K is the strike price
  'T is the maturity
  Dim ST, n As Integer, sum As Double, j As Integer
  n = UBound(S0) - LBound(S0) + 1
  Randomize
  
  For j = 1 To 5000
    'Get a series of simulated stock prices
    ST = simCorrGBMStock(S0, delta, sigma, corr, r, T)
  
    'Calculate the value of the portfolio
    Dim V As Double, i As Integer
    V = 0
    For i = 0 To n - 1
      V = V + ST(i) * w(i)
    Next i
  
    'Calculate the payoff and add to the total sum
    sum = sum + WorksheetFunction.Max(V - K, 0)
    
    DoEvents
  Next j
  
  'Output the result using risk-neutral pricing formula
  European_Basket_Call_Option_MC = Exp(-r * T) * sum / 5000
  
End Function

Sub ch6_9_ex()
  Dim S0(3) As Double, delta(3) As Double, sigma(3) As Double, corr(3, 3) As Double, w(3) As Double
  
  'initialization
  S0(0) = 70
  delta(0) = 0.02
  sigma(0) = 0.25
  corr(0, 0) = 1
  w(0) = 0.25
  
  S0(1) = 40
  delta(1) = 0
  sigma(1) = 0.3
  corr(0, 1) = 0.1
  corr(1, 0) = 0.1
  corr(1, 1) = 1
  w(1) = 0.25
  
  S0(2) = 50
  delta(2) = 0.02
  sigma(2) = 0.35
  corr(0, 2) = -0.1
  corr(2, 0) = -0.1
  corr(1, 2) = 0.3
  corr(2, 1) = 0.3
  corr(2, 2) = 1
  w(2) = 0.25
  
  S0(3) = 45
  delta(3) = 0
  sigma(3) = 0.3
  corr(0, 3) = -0.3
  corr(3, 0) = -0.3
  corr(1, 3) = 0.2
  corr(3, 1) = 0.2
  corr(2, 3) = 0.6
  corr(3, 2) = 0.6
  corr(3, 3) = 1
  w(3) = 0.25
  
  'Output the price
  MsgBox European_Basket_Call_Option_MC(S0, delta, sigma, corr, w, 0.1, 50, 1)
End Sub

Function Black_Scholes_American_Call_LS(ByVal S0 As Double, ByVal T As Double, ByVal K As Double, ByVal r As Double, ByVal delta As Double, ByVal sigma As Double, nStep As Integer) As Double
  'Use the Longstaff-Schwartz Least-Square approach to approximate
  'S0 is the initial stock price
  'T is the time to maturity
  'K is the strike price
  'r is the continuous interest rate
  'delta is the continuous dividend rate
  'sigma is the stock volatility
  'nStep is the number of steps of the binomial tree
  Randomize
  
  'Determine the number of simulations
  Dim n As Integer
  n = 10000
  
  'Calculate the time for each period
  Dim h As Double
  h = T / nStep
  
  'Create a variable to generate a complete path of all the simulated stocks
  Dim path() As Double
  ReDim path(n - 1, nStep)
  
  'Loop for the number of simulations
  Dim i As Integer, j As Integer
  For i = 0 To n - 1
  
    'Generate a complete path of the stock
    path(i, 0) = S0
    For j = 1 To nStep
      path(i, j) = simGBMStock(path(i, j - 1), sigma, r, delta, h)
    Next j
    
    DoEvents
  Next i
    
  'Variables
  Dim payoff() As Double          'hold the final payoff for all simulated paths
  Dim exercisePeriod() As Integer 'hold the exercising period
  ReDim payoff(n - 1)
  ReDim exercisePeriod(n - 1)
  
  'Look back the path
  For j = nStep To 1 Step -1
    Dim V As Collection     'continuation value
    Dim VPath As Collection 'holding path info for corresponding V
    
    Set V = New Collection
    Set VPath = New Collection
    
    Dim tempPayoff As Double
    
    'For each simulated path
    For i = 0 To n - 1
    
      'Calculate the payoff if exercise at this time point
      tempPayoff = WorksheetFunction.Max(path(i, j) - K, 0)
      
      'If it is at the last period, just put the payoff to the variable
      If j = nStep Then
        payoff(i) = tempPayoff
        exercisePeriod(i) = nStep

      'Else, exercise is optimal if payoff > 0
      ElseIf tempPayoff > 0 Then
      
        'Calculate the continuation value
        V.Add payoff(i) * Exp((j - exercisePeriod(i)) * h * r)
        VPath.Add i
        
      End If
      
      DoEvents
    Next i
    
    If j <> nStep And V.count > 0 Then

      If V.count > 3 Then
      
        'Copy the continuation values and stock price
        Dim VToReg() As Double, SToReg() As Double, S2ToReg() As Double
        ReDim VToReg(V.count - 1)
        ReDim SToReg(V.count - 1)
        ReDim S2ToReg(V.count - 1)
      
        For i = 0 To V.count - 1
          VToReg(i) = V(i + 1)
          SToReg(i) = path(VPath(i + 1), j)
          S2ToReg(i) = SToReg(i) ^ 2
          DoEvents
        Next i
      
        'Do the regression
        Dim lsSol
        lsSol = LeastSquaresSolution2(VToReg, SToReg, S2ToReg)

      End If
      
      'For each simulated path with continuation value
      For i = 0 To V.count - 1
        
        'Do the forcast
        Dim Vf As Double
        If V.count <= 3 Then Vf = V(i + 1) _
        Else Vf = lsSol(1) + lsSol(2) * SToReg(i) + lsSol(3) * S2ToReg(i)
        
        'Calculate the tempPayoff again
        tempPayoff = WorksheetFunction.Max(path(VPath(i + 1), j) - K, 0)
        
        'If the forcast value < exercise value
        If Vf < tempPayoff Then
          'exercise, record the payoff and the exercise period
          payoff(VPath(i + 1)) = tempPayoff
          exercisePeriod(VPath(i + 1)) = j
        End If
        
        DoEvents
      Next i
      
    End If
    
  Next j
  
  'Finally calculate the mean of discounted payoffs
  Dim sum As Double
  For i = 0 To n - 1
    sum = sum + payoff(i) * Exp(-exercisePeriod(i) * h * r)
  Next i
  
  Black_Scholes_American_Call_LS = sum / n

End Function

Function Black_Scholes_American_Put_LS(ByVal S0 As Double, ByVal T As Double, ByVal K As Double, ByVal r As Double, ByVal delta As Double, ByVal sigma As Double, nStep As Integer) As Double
  'Use the Longstaff-Schwartz Least-Square approach to approximate
  'S0 is the initial stock price
  'T is the time to maturity
  'K is the strike price
  'r is the continuous interest rate
  'delta is the continuous dividend rate
  'sigma is the stock volatility
  'nStep is the number of steps of the binomial tree
  Randomize
  
  'Determine the number of simulations
  Dim n As Integer
  n = 10000
  
  'Calculate the time for each period
  Dim h As Double
  h = T / nStep
  
  'Create a variable to generate a complete path of all the simulated stocks
  Dim path() As Double
  ReDim path(n - 1, nStep)
  
  'Loop for the number of simulations
  Dim i As Integer, j As Integer
  For i = 0 To n - 1
  
    'Generate a complete path of the stock
    path(i, 0) = S0
    For j = 1 To nStep
      path(i, j) = simGBMStock(path(i, j - 1), sigma, r, delta, h)
    Next j
    
    DoEvents
  Next i
    
  'Variables
  Dim payoff() As Double          'hold the final payoff for all simulated paths
  Dim exercisePeriod() As Integer 'hold the exercising period
  ReDim payoff(n - 1)
  ReDim exercisePeriod(n - 1)
  
  'Look back the path
  For j = nStep To 1 Step -1
    Dim V As Collection     'continuation value
    Dim VPath As Collection 'holding path info for corresponding V
    
    Set V = New Collection
    Set VPath = New Collection
    
    Dim tempPayoff As Double
    
    'For each simulated path
    For i = 0 To n - 1
    
      'Calculate the payoff if exercise at this time point
      tempPayoff = WorksheetFunction.Max(K - path(i, j), 0)
      
      'If it is at the last period, just put the payoff to the variable
      If j = nStep Then
        payoff(i) = tempPayoff
        exercisePeriod(i) = nStep

      'Else, exercise is optimal if payoff > 0
      ElseIf tempPayoff > 0 Then
      
        'Calculate the continuation value
        V.Add payoff(i) * Exp((j - exercisePeriod(i)) * h * r)
        VPath.Add i
        
      End If
      
      DoEvents
    Next i
    
    If j <> nStep And V.count > 0 Then

      If V.count > 3 Then
      
        'Copy the continuation values and stock price
        Dim VToReg() As Double, SToReg() As Double, S2ToReg() As Double
        ReDim VToReg(V.count - 1)
        ReDim SToReg(V.count - 1)
        ReDim S2ToReg(V.count - 1)
      
        For i = 0 To V.count - 1
          VToReg(i) = V(i + 1)
          SToReg(i) = path(VPath(i + 1), j)
          S2ToReg(i) = SToReg(i) ^ 2
          DoEvents
        Next i
      
        'Do the regression
        Dim lsSol
        lsSol = LeastSquaresSolution2(VToReg, SToReg, S2ToReg)

      End If
      
      'For each simulated path with continuation value
      For i = 0 To V.count - 1
        
        'Do the forcast
        Dim Vf As Double
        If V.count <= 3 Then Vf = V(i + 1) _
        Else Vf = lsSol(1) + lsSol(2) * SToReg(i) + lsSol(3) * S2ToReg(i)
        
        'Calculate the tempPayoff again
        tempPayoff = WorksheetFunction.Max(K - path(VPath(i + 1), j), 0)
        
        'If the forcast value < exercise value
        If Vf < tempPayoff Then
          'exercise, record the payoff and the exercise period
          payoff(VPath(i + 1)) = tempPayoff
          exercisePeriod(VPath(i + 1)) = j
        End If
        
        DoEvents
      Next i
      
    End If
    
  Next j
  
  'Finally calculate the mean of discounted payoffs
  Dim sum As Double
  For i = 0 To n - 1
    sum = sum + payoff(i) * Exp(-exercisePeriod(i) * h * r)
  Next i
  
  Black_Scholes_American_Put_LS = sum / n

End Function

Public Function LeastSquaresSolution2(y, ParamArray x())
  Dim a(), nRow As Integer, nColumn As Integer, i As Integer, j As Integer
  nRow = UBound(x(0)) + 1 'number of rows of A
  nColumn = UBound(x) + 2 'number of columns of A
  ReDim a(nRow - 1, nColumn - 1)
  
  For i = 0 To nRow - 1 'First column is 1
    a(i, 0) = 1
  Next i
  
  For j = 1 To nColumn - 1 'For each column
    For i = 0 To nRow - 1 'For each row
      a(i, j) = x(j - 1)(i)
    Next i
  Next j
  
  LeastSquaresSolution2 = WorksheetFunction.Transpose( _
                           MMult( _
                             MMult( _
                               WorksheetFunction.MInverse( _
                                 MMult( _
                                   WorksheetFunction.Transpose(a), _
                                   a _
                                 ) _
                               ), _
                               WorksheetFunction.Transpose(a) _
                             ), _
                             WorksheetFunction.Transpose(y) _
                           ) _
                         )
End Function

Public Function MMult(a, b)
  Dim result(), nRowA As Long, nColumnA As Long, nColumnB As Long, i As Long, j As Long, K As Long
  nRowA = UBound(a, 1) - LBound(a, 1) + 1
  nColumnA = UBound(a, 2) - LBound(a, 2) + 1
  nColumnB = UBound(b, 2) - LBound(b, 2) + 1
  ReDim result(nRowA - 1, nColumnB - 1)
  
  For i = 0 To nRowA - 1
    For j = 0 To nColumnB - 1
      For K = 0 To nColumnA - 1
        result(i, j) = result(i, j) + a(i + LBound(a, 1), K + LBound(a, 2)) * b(K + LBound(b, 1), j + LBound(b, 2))
      Next K
    Next j
  Next i
  
  MMult = result
End Function

Sub testBlack_Scholes_American_Put_LS()
  Cells(1, 1) = "European"
  Cells(1, 2) = Black_Scholes_Put(1, 3, 1.1, 0.06, 0, 0.15)
  Cells(3, 1) = "American Sim..."
  For i = 1 To 10
    Cells(i + 3, 1) = Black_Scholes_American_Put_LS(1, 3, 1.1, 0.06, 0, 0.15, 3)
  Next i
End Sub