devbox
Cpp五条
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

mpf4_定价欧式美式障碍Options_CRR_Leisen-Reimer_Greeks_二叉树三叉树网格_Finite differences(显式隐式)Crank-Nicolson_Imp波动率_a stock price is currently $50. over each of the n

作者：Cpp五条 | 2024-03-16 17:20:51

赞

踩

a stock price is currently $50. over each of the next two three-month period

     A derivative[dɪˈrɪvətɪv](金融)衍生工具(产品)  is a contract whose payoff depends on the value of some underlying asset. In cases where closed-form derivative pricing may be complex or even impossible, numerical procedures excel. A numerical procedure is the use of iterative computational methods in attempting to converge to a solution. One such basic implementation is a binomial tree. In a binomial tree, a node represents the state of an asset at a certain point of time associated with a price. Each node leads to two other nodes in the next time step. Similarly, in a trinomial tree, each node leads to three other nodes in the next time step. However, as the number of nodes or the time steps of trees increase, so do the computational resources that are consumed. Lattice pricing attempts to solve this problem by storing only the new information at each time step, while reusing values where possible.

     In finite difference pricing限差分定价, the nodes of the tree can also be represented as a grid. The terminal values on the grid consist of terminal conditions, while the edges of the grid represent boundary conditions in asset pricing. We will discuss the explicit method, implicit method, and the Crank-Nicolson method of the finite difference schemes to determine the price of an asset.

     Although vanilla[vəˈnɪlə] options普通期权 and certain exotics某些特殊产品 such as European barrier options欧式障碍期权 and lookback options回溯期权 can be found to have a closed-form solution, other exotic products such as Asian options do not contain a closed-form solution. In these cases, the pricing of options can be used with numerical procedures.

In this chapter, we will cover the following topics:

  • Pricing European and American options using a binomial tree
  • Using a Cox-Ross-Rubinstein binomial tree
  • Pricing options using a Leisen-Reimer tree
  • Pricing options using a trinomial tree
  • Deriving Greeks from a tree for free
  • Pricing options using a binomial and trinomial lattice
  • Finite differences with the explicit, implicit, and Crank-Nicolson method
  • Implied volatility modeling using an LR tree and the bisection method二分法

Introduction to options

     An option is a derivative of an asset that gives an owner the right but not the obligation to transact the underlying asset at a certain date for a certain price, known as the maturity date and strike price, respectively.

     A call option看涨期权 gives the buyer the right to buy an asset by a certain date for a certain price. A seller or writer of a call option is obligated to sell the underlying security to the buyer at the agreed price, should the buyer exercise his/her rights on the agreed date.

     A put option看跌期权 gives the buyer the right to sell the underlying asset by a certain date for a certain price. A seller or writer of a put option is obligated to buy the underlying security from the buyer at the agreed price, should the buyer exercise his/her rights on the agreed date.

     The most common options that are available are the European options and American options. Other exotic options include Bermudan options百慕大期权 and Asian options. This chapter will deal mainly with European and American options. A European option can only be exercised on the maturity date. An American option, on the other hand, may be exercised at any time throughout the lifetime of the option.

Binomial trees in option pricing

     In the binomial option pricing model, the underlying security at one time period, represented as a node with a given price, is assumed to traverse to two other nodes in the next time step, representing an up state and a down state. Since options are derivatives of the underlying asset, the binomial pricing model tracks the underlying conditions on a discrete-time basis. Binomial option pricing can be used to value European options, American options, as well as Bermudan options.

     The initial value of the root node is the spot price现货价格 of the underlying security with a risk-neutral probability of increase q风险中性增加概率为 q, and a risk-neutral probability of loss 1-q风险中性损失概率为 1-q, at the next time step. Based on these probabilities, the expected values of the security are calculated for each state of price increase or decrease for every time step. The terminal nodes represent every value of the expected security prices for every combination of up states and down states. We can then calculate the value of the option at every node, traverse the tree by risk-neutral expectations, and after discounting from the forward interest rates从远期利率折现后, we can derive the value of the call or put option.

Pricing European options

     Consider a two-step binomial tree. A non-dividend paying stock price starts at $50, and, in each of the two time steps, the stock may go up by 20 percent or go down by 20 percent. Suppose that the risk-free rate is five percent per annum[ˈænəm]年 and that the time to maturity, T, is two years. We would like to find the value of a European put option with a strike price K of $52. The following diagram shows the pricing of the stock and the payoffs收益 at the terminal nodes using a binomial tree:

Here, the nodes are calculated as follows:

  • Spot Price, =50
  • : The stock price today.
    p: The probability of a price rise.

    u: The factor by which the price rises (assuming it rises)

    d: The factor by which the price falls (assuming it falls).
  • Probability of price at up state : p
    the stock may go up by 20 percent : u=1.2
  • Probability of price at down state: p
    go down by 20 percent : d=0.8

  • if S (such as ) is the current price, then in the next period the price will either be  or .

At the terminal nodes, the payoff from exercising a European call option is given as follows: 

In the case of a European put option, the payoff is as follows: 

     European call and put options are usually denoted by lowercase letters, c and p, while American call and put options are usually denoted by uppercase letters, C and P.

     From the option payoff values, we can then traverse the binomial tree backward to the current time, and after discounting from the risk-free rate we will obtain our present value of the option. Traversing the tree backward takes into account the risk-neutral probabilities of the option's up states and down states.

     We may assume that investors are indifferent[ɪnˈdɪfrənt]不关心的 to risk and that expected returns on all assets are equal. In the case of investing in stocks by risk-neutral probability, the payoff from holding the stock and taking into account the up and down state possibilities would be equal to the continuously compounded risk-free rate expected in the next time step持有股票并考虑上下状态可能性的收益 将等于 下一时间步 预期的连续复合无风险利率, as follows:OR which ensures that over the small period of time Δt the expected return of the binomial model matches the expected return in a risk-neutral world, and the equation, which ensures that the variance matches.

The risk-neutral probability q of investing in the stock can be rewritten as follows: 

Are these formulas relevant to stocks? What about futures期货? 

     Unlike investing in stocks, investors do not have to make an upfront payment to take a position in a futures contract. In a risk-neutral sense, the expected growth rate from holding a futures contract is zero, and the risk-neutral probability q of investing in futures can be rewritten as follows: 

Let's calculate the risk-neutral probability q of the stock given in the preceding example:

     Consider a two-step binomial tree. A non-dividend paying stock price starts at $50, and, in each of the two time steps, the stock may go up by 20 percent or go down by 20 percent. Suppose that the risk-free rate is 5% per annum[ˈænəm]年 and that the time to maturity, T, is two years. We would like to find the value of a European put option with a strike price K of $52.

The risk-neutral probability q of investing in the stock:OR

  1. import math
  2.  
  3. r=0.05  # risk-free rate
  4. T = 2   # the time to maturity, T, is two years
  5. t = T/2 # the next time step
  6. u = 1.2 # go up by 20 percent
  7. d = 0.8 # go down by 20 percent
  8.  
  9. q = ( math.exp(r*t) - d )/(u-d)
  10. print('the risk-neutral probability q:',q)

if S (such as ) is the current price, then in the next period the price will either be  or .  

  1. Su = S0 * u
  2. Sd = S0 * d
  3. Suu = Su * u
  4. Sud = Su * d # or Sd*u
  5. Sdd = Sd * d
  6.  
  7. print( "Su:", Su )
  8. print( "Sd:", Sd )
  9. print( "Suu:", Suu )
  10. print( "Sud:", Sud )
  11. print( "Sdd:", Sdd )

  1. puu = max(0, K-Suu)
  2. pud = max(0, K-Sud)
  3. pdd = max(0, K-Sdd)
  4.  
  5. pu = math.exp(-r*t)*( puu*q + pud*(1-q) )
  6. pd = math.exp(-r*t)*( pud*q + pdd*(1-q) )
  7.  
  8. p0=math.exp(-r*t)*( pu*q + pd*(1-q) )
  9.  
  10. print( "puu:", puu)
  11. print( "pud:", pud)
  12. print( "pdd:", pdd)
  13. print( "pu: ", pu )
  14. print( "pd: ", pd )
  15. print( "p0: ", p0 )


     The payoffs of exercising the European put option at the terminal nodes are $0, $4, and $20 at the respective states. The present value of the put option can be priced as follows:
<==strike price K : $52

  1. pu=math.exp(-r*t)*( 0*q + 4*(1-q) )
  2. pu

vs

  1. pd=math.exp(-r*t)*( 4*q + 20*(1-q) )
  2. pd

 

  1. p0=math.exp(-r*t)*( pu*q + pd*(1-q) )
  2. p0

 

  1. pt=math.exp(-r*T)*( 0*q**2 + 4*q*(1-q)*2 + 20*(1-q)**2 )
  2. pt


<==

######################## Call option

A stock price is currently $50. Over each of the next two three-month periods it is expected to go up by 6% or down by 5%. The risk-free interest rate is 5% per annum with continuous compounding. What is the value of a six-month European call option with a strike price of $51?Figure S12.1
     A tree describing the behavior of the stock price is shown in Figure S12.1. The risk-neutral probability of an up move, p, is given by

  1. S0 = 50  # A stock price is currently $50
  2. K = 51   # a strike price of $51
  3. r=0.05   # risk-free rate
  4. T = 6/12 # the time to maturity, T, in years
  5. n_t = 2 # the number of time steps
  6. t = T/n_t
  7. u = 1.06 # it is expected to go up by 6% 
  8. d = 0.95 # it is expected to go down by 5% 
  9.  
  10. q = ( math.exp(r*t) - d )/(u-d)
  11. print('the risk-neutral probability q:',q)

 

if S (such as ) is the current price, then in the next period the price will either be  or . 

  1. Su = S0 * u
  2. Sd = S0 * d
  3. Suu = Su * u
  4. Sud = Su * d # or Sd*u
  5. Sdd = Sd * d
  6.  
  7. print( "Su:", Su )
  8. print( "Sd:", Sd )
  9. print( "Suu:", Suu )
  10. print( "Sud:", Sud )
  11. print( "Sdd:", Sdd )

  1. puu = max(0, Suu-K)
  2. pud = max(0, Sud-K)
  3. pdd = max(0, Sdd-K)
  4.  
  5. pu = math.exp(-r*t)*( puu*q + pud*(1-q) )
  6. pd = math.exp(-r*t)*( pud*q + pdd*(1-q) )
  7.  
  8. p0=math.exp(-r*t)*( pu*q + pd*(1-q) )
  9.  
  10. print( "puu:", puu)
  11. print( "pud:", pud)
  12. print( "pdd:", pdd)
  13. print( "pu: ", pu )
  14. print( "pd: ", pd )
  15. print( "p0: ", p0 )


     There is a payoff from the call option of 56.18-51=5.18 for the highest final node (which corresponds to two up moves) zero in all other cases. The value of the option is therefore

  1. pt=math.exp(-r*T)*( puu*q**2 + pud*q*(1-q)*2 + pdd*(1-q)**2 )
  2. pt



######################## 

Writing the StockOption base class 

     Before going any further and implementing the various pricing models that we are about to discuss, let's create a StockOption class to store and calculate the common attributes of the stock option that will be reused throughout this chapter: 

  1. import math
  2. # Stores common attributes of stock option
  3. class StockOption(object):
  4.     def __init__(
  5.         self, S0, K, r =0.05, T=1, N=2, pu=0, pd=0,
  6.         div=0, sigma=0, is_put = False, is_am=False
  7.     ):
  8.         """
  9.         Initialize the stock option base class.
  10.         Defaults to European call unless specified.
  11.         : param S0: initial stock price
  12.         : param K: strike price
  13.         : param r: risk-free interest rate
  14.         : param T: time to maturity
  15.         : param N: number of time steps
  16.         : param pu: probability at up state
  17.         : param pd: probability at down state
  18.         : param div: Dividend yield
  19.         : param is_put: True for a put option,
  20.                         False for a call option
  21.         : param is_am: True for an American option,
  22.                        False for a European option
  23.         """
  24.         self.S0=S0
  25.         self.K=K
  26.         self.r=r
  27.         self.T=T
  28.         self.N = max(1,N)
  29.         self.STs = [] # Declare the stock prices tree
  30.         
  31.         # Optional parameters used by derived classes
  32.         self.pu, self.pd = pu, pd
  33.         self.div = div
  34.         self.sigma = sigma
  35.         self.is_call = not is_put
  36.         self.is_european = not is_am
  37.         # @property:https://www.programiz.com/python-programming/property
  38.     @property
  39.     def dt(self):
  40.         # Single time step, in years
  41.         return self.T/float(self.N)
  42.     
  43.     @property
  44.     def df(self):
  45.         # The discount factor
  46.         return math.exp( -(self.r - self.div)*self.dt )

  1. s = StockOption(50, 52,T=2)
  2. s.df #e^(-rt)

 == code: math.exp(-r*1) and r=0.05

     The

  • current underlying price(S0),
  • strike price(K),
  • risk-free rate(r),
  • time to maturity(T),
  • and number of time steps(N)

are compulsory common attributes for pricing options. The

  • self.div = 0
  • delta of the time step dt (=self.T/float(self.N) in years)
  • and the discount factor df ( =math.exp( -(self.r - self.div)*self.dt ) ) ==>

are computed as properties of the class and may be overwritten by implementing classes if needed.

A class for European options using a binomial tree

     The Python implementation of the binomial option pricing model of a European option is given as the BinomialEuropeanOption class, inheriting the common attributes of the option from the StockOption class. The implementations of the methods in this class are as follows:

  • 1.The price() method of the BinomialEuropeanOption class is the entry point for all the instances of this class
  • 2. It invokes the setup_parameters() method to set up the required model parameters, and then calls the init_stock_price_tree() method to simulate the expected values of the stock prices for the period up until maturity (T)
  • 3. Finally, the begin_tree_traversal() method is called to initialize the payoff array and store the discounted payoff values, as it traverses the binomial tree back to the present time
  • 4. The payoff tree nodes are returned as a NumPy array object, where the present value of the European option is found at the initial node

The class implementation of BinomialEuropeanOption is given in the following Python code:

  1. import math
  2. import numpy as np
  3. from decimal import Decimal
  4.  
  5. # Price a European option by the binomial tree model
  6.  
  7. class BinomialEuropeanOption(StockOption):
  8.     def setup_parameters(self):
  9.         # Required calculations for the model
  10.         self.M = self.N+1 # Number of terminal nodes of tree
  11.         
  12.         self.u = 1 + self.pu # Expected value in the up state, e.g 1+0.2=1.2
  13.         self.d = 1 - self.pd # Expected value in the down state e.g 1-0.2=0.8
  14.         
  15.         # qu : The risk-neutral (up state) probability 
  16.         self.qu = ( math.exp( (self.r-self.div)*self.dt )-self.d 
  17.                   )/(self.u-self.d)
  18.         # qd : The risk-neutral (down state) probability
  19.         self.qd = 1-self.qu
  20.         
  21.     def init_stock_price_tree( self ):
  22.         # Initialize terminal price nodes to zeros
  23.         self.STs = np.zeros( self.M )
  24.         
  25.         # Calculate expected stock prices for each terminal node
  26.         for i in range( self.M ):
  27.             self.STs[i] = self.S0 * ( self.u**(self.N-i) )*( self.d**i )
  28.         # Suu = $72, Sud=Sdu=$48, Sdd=$32
  29.     
  30.     def init_payoffs_tree( self ):
  31.         # Returns the payoffs when the option expires at terminal nodes
  32.         if self.is_call: # call option 
  33.             return np.maximum( 0, self.STs - self.K )
  34.         else: # put option 
  35.             return np.maximum( 0, self.K - self.STs )
  36.             # Puu=$0, Pud=Pdu=$4, Pdd=$20 
  37.     
  38.     def traverse_tree( self, payoffs ):
  39.         """
  40.         Starting from the time the option expires, traverse
  41.         backwards and calculate discounted payoffs at each node
  42.         """
  43.         # N: number of time steps
  44.         # payoffs : [Puu=$0, Pud=Pdu=$4, Pdd=$20]
  45.         for i in range( self.N ):
  46.             payoffs = ( payoffs[:-1] * self.qu + 
  47.                         payoffs[1:] * self.qd
  48.                       ) * self.df # df: The discount factor
  49.                       # * math.exp( -(self.r - self.div)*self.dt )
  50.             # array([1.41475309, 9.46393007])
  51.             # array([4.19265428])
  52.         return payoffs
  53.     
  54.     def begin_tree_traversal( self ):
  55.         # Initialize the payoffs at the terminal node
  56.         payoffs = self.init_payoffs_tree() # [Puu=$0, Pud=Pdu=$4, Pdd=$20] 
  57.         return self.traverse_tree( payoffs )
  58.     
  59.     ###### Entry point of the pricing implementation ######   
  60.     def price( self ):
  61.         self.setup_parameters()
  62.         self.init_stock_price_tree() # STs : Suu = $72, Sud=Sdu=$48, Sdd=$32
  63.         payoffs = self.begin_tree_traversal()
  64.         # Option value converges to first node
  65.         return payoffs[0]

     Let's take the values from the two-step binomial tree example we discussed earlier to price the European put option:

  1. eu_option = BinomialEuropeanOption(
  2.     50, 52, r=0.05, T=2, N=2, pu=0.2, pd=0.2, is_put=True)
  3. print( 'European put option price is:', eu_option.price() )


     Using the binomial option pricing model gives us a present value of $4.19 for the European put option.

  1. eu_option = BinomialEuropeanOption(
  2.     50, 51, r=0.05, T=6/12, N=2, pu=0.06, pd=0.05, is_put=True)
  3. print( 'European put option price is:', eu_option.price() )


###############

Call option

  1. eu_option = BinomialEuropeanOption(
  2.     50, 51, r=0.05, T=6/12, N=2, pu=0.06, pd=0.05, is_put=False)
  3.  
  4. print( 'European Call option price is:', eu_option.price() )
  5. print( 'Expected stock prices for each terminal node:', eu_option.STs )
  6. print( "Expected payoff for each terminal node:",
  7.        eu_option.init_payoffs_tree() 
  8.      )

 

put option

  1. eu_option = BinomialEuropeanOption(
  2.     50, 51, r=0.05, T=6/12, N=2, pu=0.06, pd=0.05, is_put=True)
  3.  
  4. print( 'European put option price is:', eu_option.price() )
  5. print( 'Expected stock prices for each terminal node:', eu_option.STs )
  6. print( "Expected payoff for each terminal node:",
  7.        eu_option.init_payoffs_tree() 
  8.      )

 

Verify that the European call and European put prices satisfy put–call parity看跌期权平价?

 The value of the put plus the stock price is
 The value of the call plus the present value of the strike price is

eu_option.price() + eu_option.K * math.exp(-eu_option.r * eu_option.T )

=  This verifies that put–call parity holds

     To test whether it worth exercising the option early we compare the value calculated for the option at each node with the payoff from immediate exercise. At node C the payoff from immediate exercise is  . Because this is greater than 2.8664, the option should be exercised at this node. The option should not be exercised at either node A or node B.  

 ###############

A class for American options using a binomial tree

     Unlike European options, which can only be exercised at maturity, American options can be exercised at any time during their lifetime.

     To implement the pricing of American options in Python, in the same way we did with the BinomialEuropeanOption class, create a class named BinomialTreeOption that inherits the Stockoption class. The parameters that are used in the setup_parameters() method remain the same, except for the removal of an unused M parameter(Number of terminal nodes of tree).

The methods that are used in American options are as follows:  K=52

  • init_stock_price_tree : Uses a two-dimensional NumPy array to store the expected returns of the stock prices for all time steps. This information is used to calculate the payoff values from exercising the option at each period. This method is written as follows:

    # self.STs : [array([50]), array([60., 40.]), array([72., 48., 32.])]

    if S (such as ) is the current price, then in the next period the price will either be  or .

    The up and down factors are calculated using the underlying volatility, , and the time duration of a step, t, measured in years (using the day count convention of the underlying instrument). From the condition that the variance of the log of the price is , we have:, but here we give u =1.2, d=0.8
  • init_payoffs_tree : Creates the payoff tree as a two-dimensional NumPy array, starting with the intrinsic values of the option at maturity. This method is written as follows:
  • check_early_exercise : Returns the maximum payoff values between exercising the American option early提前行使美式期权 and not exercising the option at all完全不行使期权. This method is written as follows:
  • traverse_tree : e^(−rt) × (q×P_up​+(1−q)×P_down​) ==>[1.41475309 9.46393007] ==>This also includes the invocation of the check_early_exercise() method to check whether it is optimal to exercise an American option early at every time step. ==>np.maximum( payoffs, self.STs[node]-self.K ) or np.maximum( payoffs, self.K-self.STs[node] ) ==>[1.41475309 12.]
    [1.41475309     12]<==np.maximum( payoffs, self.K-self.STs[node] )<==

    The following formula to compute the expectation value is applied at each node:==>
  • The implementation of the begin_tree_traversal() ==>payoffs: [ 0. 4. 20.]
    and the price() methods remains the same.
  • The BinomialTreeOption class can price both European and American options when the is_put keyword argument is set to False or True during instantiation of the class, respectively.

 The following code is for pricing the American option:

  1. import math
  2. import numpy as np
  3.  
  4. # Price a European or American option by the binomial tree
  5. class BinomialTreeOption(StockOption):
  6.     def setup_parameters( self ):
  7.         self.u = 1 + self.pu # Expected value % in the up state, e.g 1+0.2=1.2
  8.         self.d = 1 - self.pd # Expected value % in the down state e.g 1-0.2=0.8
  9.         # qu : The risk-neutral (up state) probability 
  10.         self.qu = ( math.exp( (self.r-self.div)*self.dt )-self.d 
  11.                   )/(self.u-self.d)
  12.         # qd : The risk-neutral (down state) probability
  13.         self.qd = 1-self.qu
  14.         
  15.     def init_stock_price_tree(self):
  16.         # Initialize a 2D tree at T=0
  17.         # to store the expected returns of the stock prices for all time steps
  18.         self.STs = [ np.array([self.S0]) ] # array([50]) # root
  19.         
  20.         # Simulate the possible stock prices path
  21.         # calculate the payoff values 
  22.         # from exercising the option at each period
  23.         for i in range( self.N ): # N: number of time steps(branch_nodes)
  24.             prev_branches = self.STs[-1] # at last time step
  25.             st = np.concatenate(             # note, cross node: Sud = Sdu
  26.                 ( prev_branches*self.u,      # go up
  27.                   [prev_branches[-1]*self.d] # the last one goes down
  28.                 )# array([60., 40.])
  29.             )    # array([72., 48., 32.])
  30.             
  31.             self.STs.append(st) # Add nodes at each time step
  32.     
  33.     # the intrinsic values of the option at maturity
  34.     def init_payoffs_tree(self):
  35.         if self.is_call:# call option
  36.             return np.maximum( 0, self.STs[self.N]-self.K )
  37.         else:# put option
  38.             return np.maximum( 0, self.K-self.STs[self.N] )
  39.              # e.g. payoffs: [ 0.  4. 20.]
  40.     
  41.     # Returns the maximum payoff values between 
  42.     # exercising the American option early 
  43.     # and not exercising the option at all
  44.     def check_early_exercise( self, payoffs, node ):
  45.         if self.is_call:
  46.             return np.maximum( payoffs, self.STs[node]-self.K )
  47.         else:
  48.             # payoffs: [1.41475309 9.46393007]
  49.             # payoffs at maturity: [-8. 12.]
  50.             # np.maximum ==>[1.41475309 12.]
  51.             return np.maximum( payoffs, self.K-self.STs[node] )
  52.     
  53.     def traverse_tree( self, payoffs ):
  54.         """
  55.         Starting from the time the option expires, traverse
  56.         backwards and calculate discounted payoffs at each node
  57.         """
  58.         for i in reversed( range(self.N) ):
  59.             # The payoffs from NOT exercising the option
  60.             payoffs = ( payoffs[:-1]*self.qu +
  61.                         payoffs[1:]*self.qd
  62.                       ) * self.df # df: The discount factor
  63.                       # * math.exp( -(self.r - self.div)*self.dt )
  64.             print("At the time step",i,", payoffs:",payoffs)
  65.             # [1.41475309 9.46393007]
  66.             # [4.19265428]
  67.             # Payoffs from exercising, for American options
  68.             if not self.is_european:
  69.                 payoffs = self.check_early_exercise( payoffs, i )
  70.                 print("At the time step",i,", max payoffs:",payoffs)
  71.                 
  72.         return payoffs
  73.     
  74.     def begin_tree_traversal( self ):
  75.         # the intrinsic values of the option at maturity
  76.         payoffs = self.init_payoffs_tree() # e.g. payoffs: [ 0. 4. 20.]
  77.         print("the payoffs at maturity", payoffs)
  78.         return self.traverse_tree( payoffs ) # calculate option values
  79.     
  80.     ###### Entry point of the pricing implementation ######
  81.     def price(self):
  82.         self.setup_parameters()
  83.         self.init_stock_price_tree()
  84.         # self.STs : [array([50]), array([60., 40.]), array([72., 48., 32.])]
  85.         print("Nodes:",self.STs)
  86.         payoffs = self.begin_tree_traversal()
  87.         #print(payoffs)
  88.         return payoffs[0]

  1. am_option = BinomialTreeOption(50, 52, 
  2.     r=0.05, T=2, N=2, pu=0.2, pd=0.2, is_put=True, is_am=True)
  3.  
  4. print( "American put option price is:", am_option.price() )


     The American put option is priced at $5.0896. Since American options can be exercised at any time and European options can only be exercised at maturity, this added flexibility of American options increases their value over European options in certain circumstances.

################

European put option

  1. eu_option = BinomialTreeOption(
  2.     50, 51, r=0.05, T=6/12, N=2, pu=0.06, pd=0.05, is_put=True, is_am=False)
  3.  
  4. print( "European put option price is:", eu_option.price() )

American put option

  1. am_option = BinomialTreeOption(
  2.     50, 51, r=0.05, T=6/12, N=2, pu=0.06, pd=0.05, is_put=True, is_am=True)
  3.  
  4. print( "American put option price is:", am_option.price() )


     The American put option is priced at $1.6456. Since American options can be exercised at any time and European options can only be exercised at maturity, this added flexibility of American options increases their value over European options in certain circumstances.

     To test whether it worth exercising the option early we compare the value calculated for the option at each node with the payoff from immediate exercise. At node C the payoff from immediate exercise is . Because this is greater than 2.8664, the option should be exercised at this node. The option should not be exercised at either node A or node B.  
################

     For American call options on an underlying asset that does not pay dividends, there might not be an extra value over its European call option counterpart. Because of the time value of money, it costs more to exercise the American call option today before the expiration at the strike price than at a future time with the same strike price.今天在到期前以行使价行使美式看涨期权的成本要比在未来以相同行使价行使更高的成本。 For an in-the-money American call option, exercising the option early losses

  • the benefit of protection against adverse price movement不利价格变动 below the strike price,
  • as well as its intrinsic time value.
  • With no entitlement of dividend payments, there are no incentives to exercise American call options early.由于没有支付股息的权利，因此没有动力提前行使美式看涨期权。

  1. am_option2 = BinomialTreeOption(10, 12,
  2.     r=0.006, T=9, N=3, pu=0.2, pd=0.2, is_put=True, is_am=False)
  3.  
  4. print("European put option price is:", am_option2.price())

  1. am_option2 = BinomialTreeOption(10, 12,
  2.     r=0.006, T=9, N=3, pu=0.2, pd=0.2, is_put=True, is_am=True)
  3.  
  4. print("American put option price is:", am_option2.price())

 

The Cox–Ross–Rubinstein model

     In the preceding examples, we assumed that the underlying stock price would increase by 20 percent and decrease by 20 percent in the respective u>=1 up state and 0<d<=1 down state. The Cox-Ross-Rubinstein (CRR) model proposes that, over a short period of time in the risk-neutral world, the binomial model matches the mean and variance of the underlying stock. The volatilityof the underlying stock, or the standard deviation of returns of the stock, is taken into account as follows:http://home.cerge-ei.cz/petrz/fm/f400n10.pdf
From the condition that the variance of the log of the price is ,

######################

     Calculate u, d, and p when a binomial tree is constructed to value an option on a foreign currency. The tree step size is one month(delta t =1/12), the domestic interest rate is 5% per annum, the foreign interest rate is 8% per annum, and the volatility is 12% per annum.==>

 and q=0
==>

     The volatility of a non-dividend-paying stock whose price is $78, is 30%. The risk-free

rate is 3% per annum (continuously compounded) for all maturities. Calculate values for

u, d, and p when a two-month time step is used(2/12=).

==>sigma=0.03
q=0 ==>==>
######################
==>

A class for the CRR binomial tree option pricing model

     The implementation of the binomial CRR model remains the same as the binomial tree we discussed earlier, with the exception of the u and d model parameters.

     In Python, let's create a class named BinomialCRROption and simply inherit the BinomialTreeOption class. Then, all that we need to do is override the setup_parameters() method with values from the CRR model.

     Instances of the BinomialCRROption object will invoke the price() method, which invokes all other methods of the parent BinomialTreeOption class, except the overwritten setup_parameters() method:
and

  1. import math
  2.  
  3. # Price an option by the binomial CRR model
  4.  
  5. class BinomialCRROption( BinomialTreeOption ):
  6.     def setup_parameters(self):  # self.dt = self.T/float(self.N)
  7.         # Expected stock price % in the up state  
  8.         self.u = math.exp( self.sigma * math.sqrt(self.dt) )
  9.         # Expected stock price % in the down state
  10.         self.d = 1./self.u
  11.  
  12.         # Calculate the transition probabilities
  13.         # qu : The risk-neutral (up state) probability 
  14.         self.qu = ( math.exp( (self.r-self.div)*self.dt )-self.d 
  15.                   )/(self.u-self.d)
  16.         # qd : The risk-neutral (down state) probability
  17.         self.qd = 1-self.qu

     Again, consider the two-step binomial tree. The non-dividend paying stock has a current price of $50 and a volatility of 30 percent. Suppose that the risk-free rate r is five percent per annum and the time to maturity T is two years. We would like to find the value of a European put option with a strike price K of $52 by the CRR model:

  1. eu_option = BinomialCRROption( 50, 52, 
  2.                 r=0.05, T=2, N=2, sigma=0.3, is_put=True )
  3. print('European put option price is:', eu_option.price() )

 

  1. am_option = BinomialCRROption( 50, 52,
  2.                 r=0.05, T=2, N=2, sigma=0.3, is_put=True, is_am=True)
  3. print( "American put option price is:", am_option.price() )

 

     By using the CRR two-step binomial tree model, the price of the European put option and the American put option are $6.2457 and $7.4284, respectively.

######################

     The volatility of a non-dividend-paying stock whose price is $78, is 30%. The risk-free

rate is 3% per annum (continuously compounded) for all maturities. Calculate values for

u, d, and p when a two-month time step is used(2/12=).

What is the value of a four-month European call option with a strike price of $80 given by a two-step binomial tree

==>
q=0 ==>==>

T=2/12*2

  1. eu_option = BinomialCRROption( 78, 80, 
  2.                 r=0.03, T=2/12*2, N=2, sigma=0.3, is_put=False )
  3. print('European put option price is:', eu_option.price() )

     Suppose a trader sells 1,000 options (10 contracts). What position in the stock is necessary to hedge the trader’s position at the time of the trade?

     The value of the option is $4.67. The initial delta is 9.58/(88.16 – 69.01) which is almost exactly 0.5 so that 500 shares should be purchased.

Problem 12.19

The current price of a non-dividend-paying biotech stock is $140 with a volatility of 25%. The risk-free rate is 4%. For a three-month time step:

  • (a) What is the percentage up movement?
  • (b) What is the percentage down movement?
  • (c) What is the probability of an up movement in a risk-neutral world?
  • (d) What is the probability of a down movement in a risk-neutral world?

Use a two-step tree to value a six-month European call option and a six-month European put option. In both cases the strike price is $150.

  1. == 1.1331. The percentage up movement is 13.31%
  2. d = 1/u = 0.8825. The percentage down movement is 1-d = 11.75%
  3. The probability of an up movement is
  4. The probability of a down movement is 0.4911=1-p.

The tree for valuing the call is in Figure S12.6a and that for valuing the put is in Figure S12.6b. The values are 7.56 and 14.58, respectively. 

  1. eu_call_option = BinomialCRRLattice(140, 150, r=0.04, T=6/12, 
  2.                       sigma=0.25, N=2, is_put=False)
  3. print("call option price:",eu_call_option.price())

 

  1. eu_put_option = BinomialCRRLattice(140, 150, r=0.04, T=6/12, 
  2.                       sigma=0.25, N=2, is_put=True)
  3. print("put option price:",eu_put_option.price())

 

######################

Using a Leisen-Reimer tree

     In the binomial models we discussed earlier, we made several assumptions about the probability of up and down states, as well as the resulting risk-neutral probabilities. Besides the binomial model with CRR(Cox-Ross-Rubinstein) parameters that we discussed, other forms of parameterization that are discussed widely in mathematical finance include the Jarrow-Rudd parameterization, Tian parameterization, and Leisen-Reimer parameterization. Let's take a look at the Leisen-Reimer model in detail.

     Dr. Dietmar Leisen and Matthias Reimer proposed a binomial tree model with the purpose of approximating to the Black-Scholes solution(https://blog.csdn.net/Linli522362242/article/details/96593981) as the number of step increases. It is known as the Leisen-Reimer (LR) tree, and the nodes do not recombine at every alternate step. It uses an inversion formula to achieve better accuracy during tree traversal.节点不会在每个交替步骤中重新组合。 它使用反演公式在树遍历期间实现更好的准确性。

     A detailed explanation of the formula is given in the paper Binomial Models For Option Valuation - Examining And Improving Convergence, March 1995, which is available at https://papers.ssrn.com/sol3/papers.cfm?abstract_id=5976. We will be using method two of the Peizer and Pratt inversion function f with the following characteristic parameters:  
             



     The

  • parameter is the current stock price,
  • K is the strike price of the option,
  • σ is the annualized volatility of the underlying stock, (volatility in annualized termshttps://blog.csdn.net/Linli522362242/article/details/122955700)
  • T is the time to maturity of the option,
  • r is the annualized risk-free rate,
  • y is the dividend yield,
  • and Δt is the time interval between each tree step.

A class for the LR binomial tree option pricing model

     The Python implementation of the LR tree is given in the following BinomialLROption class. Similar to the BinomialCRROption class, we simply inherit the BinomialTreeOption class and override the variables in the setup_parameters method with those of the LR tree model: 

  1. import math
  2.  
  3. # Price an option by the Leisen-Reimer tree
  4.  
  5. class BinomialLROption( BinomialTreeOption ):
  6.     def pp_2_inversion( self, z, n ):
  7.         # copysign : return the value of the first parameter and the sign of the second parameter
  8.         # math.copysign(4, -1)     ==> -4.0
  9.         # math.copysign(-8, 97.21) ==> 8.0
  10.         # math.copysign(-43, -76)  ==> -43.0
  11.         return 0.5 + math.copysign(1,z) *\
  12.             math.sqrt( 0.25-0.25*math.exp( 
  13.                         -( ( z/( n+1./3.+0.1/(n+1) ) )**2.
  14.                          ) * ( n+1./6. ) )
  15.                      )
  16.     
  17.     def setup_parameters( self ):
  18.         odd_N = self.N if (self.N%2==0) else (self.N+1)
  19.         d1 = ( math.log(self.S0/self.K) + 
  20.               ( (self.r-self.div) + (self.sigma**2)/2. )*self.T
  21.              ) / ( self.sigma * math.sqrt(self.T) )
  22.  
  23. #         d2 = ( math.log(self.S0/self.K) + 
  24. #               ( (self.r-self.div) - (self.sigma**2)/2. )*self.T
  25. #              ) / ( self.sigma * math.sqrt(self.T) )
  26.         d2 = d1-self.sigma*math.sqrt(self.T)
  27.         
  28.         pbar = self.pp_2_inversion( d1, odd_N )
  29.         self.p = self.pp_2_inversion( d2, odd_N )
  30.         
  31.         # self.dt : self.T/float(self.N)
  32.         # self.df : math.exp( -(self.r - self.div)*self.dt )
  33.         # self.u : Expected value in the up state, e.g 1+0.2=1.2
  34.         self.u = 1/self.df * pbar/self.p
  35.         # self.d : Expected value in the down state e.g 1-0.2=0.8
  36.         self.d = (1/self.df - self.p*self.u)/(1-self.p)
  37.         
  38.         # qu : The risk-neutral (up state) probability
  39.         self.qu = self.p
  40.         # qd : The risk-neutral (down state) probability
  41.         self.qd = 1-self.p

Using the same examples that we used previously, we can price the options using an LR tree: 

  1. eu_option = BinomialLROption(
  2.     50, 52, r=0.05, T=2, N=2, sigma=0.3, is_put=True)
  3.  
  4. print('European put option price is:', eu_option.price() )

  1. eu_option = BinomialLROption(
  2.     50, 52, r=0.05, T=2, N=4, sigma=0.3, is_put=True)
  3.  
  4. print('European put option price is:', eu_option.price() )

 

  1. am_option = BinomialLROption(
  2.     50, 52, r=0.05, T=2, N=4, sigma=0.3, is_put=True, is_am=True )
  3.  
  4. print('American put option price is:', am_option.price() )

 
     By using the LR binomial tree model with four time steps, the price of the European put option and the American put option are $5.87865 and $6.7636, respectively.

The Greeks for free

     In the binomial tree pricing models that we have covered so far, we traversed up and down the tree at each point in time to determine the node values. From the information at each node, we can reuse these computed values easily. One such use is the computation of Greeks. 

     The Greeks measure the sensitivities of the price of derivatives, such as options with respect to changes in the parameters of the underlying asset, often represented by Greek letters. In mathematical finance, the common names associated with Greeks include alpha, beta, delta, gamma, vega, theta, and rho.

     Two particularly useful Greeks for options are delta and gamma.

  • Delta measures the sensitivity of the option price with respect to the underlying asset price.
  • Gamma measures the rate of change in delta with respect to the underlying price.

     As shown in the following diagram, an additional layer of nodes is added around our original 2-step tree to make it a 4-step tree, which extends 2 steps backward in time. Even with additional terminal payoff nodes, all nodes will contain the same information as our original 2-step tree. Our option value of interest is now located in the middle of the tree at t=0:
==>
     Notice that at t=0 there exists two additional nodes' worth of information that we can use to compute the delta formula, as follows:

  • : option value in the up state at time t=0
  • :option value in the up state at time t=0
  • option value: 
  •  or :respective stock price in the up state at time t=0
  •  or :respective stock price in the down state at time t=0 

     The delta formula states that the difference in the option prices in the up and down state is represented as a unit of the difference between the respective stock prices at time t=0.

Conversely, the gamma formula can be computed as follows:

     The gamma formula states that the difference of deltas between the option prices in the up node and the down node against the initial node value are computed as a unit of the differences in price of the stock at the respective states. 

A class for Greeks with the LR binomial tree

     To illustrate the computation of Greeks with the LR tree(Leisen-Reimer tree), let's create a new class named BinomialLRWithGreeks that inherits the BinomialLROption class with our own implementation of the price method. 

     In the price method, we will start by calling the setup_parameters() method of the parent class to initialize all variables required by the LR tree. However, this time, we will also call the new_stock_price_tree() method, which is a new method that's used to create an extra layer of nodes around the original tree.

     The begin_tree_traversal() method is called to perform the usual LR tree implementation in the parent class. The returned NumPy array object now contains information on the three nodes at t=0, where the middle node is the option price. The payoffs in the up and down states at t=0 are in the first and last index of the array, respectively.

     With this information, the price() method computes and returns the option price, the delta, and the gamma values together:

  1. import numpy as np
  2.  
  3. # Compute option price, delta and gamma by the LR tree
  4.  
  5. class BinomialLRWithGreeks( BinomialLROption ):
  6.     def new_stock_price_tree( self ):
  7.         """
  8.         Creates an additional layer of nodes to our
  9.         original stock price tree
  10.         """
  11.         self.STs = [np.array([self.S0 * self.u/self.d,
  12.                               self.S0,
  13.                               self.S0 * self.d/self.u
  14.                              ]),# at t=0
  15.                    ]
  16.         
  17.         # Simulate the possible stock prices path
  18.         # calculate the payoff values 
  19.         # from exercising the option at each period   
  20.         for i in range( self.N ): # N : number of time steps
  21.             prev_branches = self.STs[-1]
  22.                                   # note, cross node: Sud = Sdu
  23.             st = np.concatenate( (prev_branches*self.u,     # branch_nodes : go up
  24.                                   [prev_branches[-1]*self.d]# the last one goes down
  25.                                  ) )
  26.             self.STs.append(st)   # Add nodes at each time step
  27.             
  28.     def price(self):
  29. #      def setup_parameters( self ):
  30. #         odd_N = self.N if (self.N%2==0) else (self.N+1)
  31. #         d1 = ( math.log(self.S0/self.K) + 
  32. #               ( (self.r-self.div) + (self.sigma**2)/2. )*self.T
  33. #              ) / ( self.sigma * math.sqrt(self.T) )
  34. #         d2 = d1-self.sigma*math.sqrt(self.T)
  35.         
  36. #         pbar = self.pp_2_inversion( d1, odd_N )
  37. #         self.p = self.pp_2_inversion( d2, odd_N )
  38.         
  39. #         # self.dt : self.T/float(self.N)
  40. #         # self.df : math.exp( -(self.r - self.div)*self.dt )
  41. #         # self.u : Expected value in the up state, e.g 1+0.2=1.2
  42. #         self.u = 1/self.df * pbar/self.p
  43. #         # self.d : Expected value in the down state e.g 1-0.2=0.8
  44. #         self.d = (1/self.df - self.p*self.u)/(1-self.p)
  45.         
  46. #         # qu : The risk-neutral (up state) probability
  47. #         self.qu = self.p
  48. #         # qd : The risk-neutral (down state) probability
  49. #         self.qd = 1-self.p       
  50.         self.setup_parameters()
  51.         self.new_stock_price_tree()
  52.         print("Nodes:",self.STs)
  53.         
  54. #     # the intrinsic values of the option at maturity
  55. #     def init_payoffs_tree(self):
  56. #         if self.is_call:# call option
  57. #             return np.maximum( 0, self.STs[self.N]-self.K )
  58. #         else:# put option
  59. #             return np.maximum( 0, self.K-self.STs[self.N] )
  60. #              # e.g. payoffs: [ 0.  4. 20.]
  61.  
  62. #     def traverse_tree( self, payoffs ):
  63. #         """
  64. #         Starting from the time the option expires, traverse
  65. #         backwards and calculate discounted payoffs at each node
  66. #         """
  67. #         for i in reversed( range(self.N) ):
  68. #             # The payoffs from NOT exercising the option
  69. #             payoffs = ( payoffs[:-1]*self.qu +
  70. #                         payoffs[1:]*self.qd
  71. #                       ) * self.df # df: The discount factor
  72. #                       # * math.exp( -(self.r - self.div)*self.dt )
  73. #             print("At the time step",i,", payoffs:",payoffs)
  74. #             # [1.41475309 9.46393007]
  75. #             # [4.19265428]
  76. #             # Payoffs from exercising, for American options
  77. #             if not self.is_european:
  78. #                 payoffs = self.check_early_exercise( payoffs, i )
  79. #                 print("At the time step",i,", max payoffs:",payoffs)                
  80. #         return payoffs
  81.  
  82. #     def begin_tree_traversal( self ):
  83. #         # the intrinsic values of the option at maturity
  84. #         payoffs = self.init_payoffs_tree() # e.g. payoffs: [ 0. 4. 20.]
  85. #         print("the payoffs at maturity", payoffs)
  86. #         return self.traverse_tree( payoffs )
  87.  
  88.         payoffs = self.begin_tree_traversal()
  89.     
  90.         #option value is now in the middel node at t=0
  91.         # self.STs[0] : [ S0*u/d , S0, S0*d/u ]
  92.         option_value = payoffs[ len(payoffs)//2 ] #option_value at S0(stock initial price)
  93.         
  94.         payoff_up = payoffs[0]    # at S0*u/d : option value in the up state
  95.         payoff_down = payoffs[-1] # at S0*d/u : option value in the down state
  96.         
  97.         S_up = self.STs[0][0]     # S0*u/d : stock price in the up state
  98.         S_down = self.STs[0][-1]  # S0*d/u : stock price in the down state
  99.         
  100.         dS_up = S_up - self.S0      # stock price in the up state - S0 
  101.         dS_down = self.S0 - S_down  # S0 - stock price in the down state
  102.         
  103.         # Calculate the delta value
  104.         # the difference in the 'option prices' in the up and down state 
  105.         # is represented as a unit of 
  106.         # the difference between the respective 'stock prices' at time t=0.
  107.         dV = payoff_up - payoff_down
  108.         dS = S_up - S_down
  109.         delta = dV/dS
  110.         
  111.         # Calculate the gamma value
  112.         gamma = ( (payoff_up - option_value)/dS_up
  113.                   - ( option_value - payoff_down) /dS_down
  114.                 )/( (self.S0 + S_up)/2 - (self.S0 + S_down)/2 )
  115.         
  116.         return option_value, delta, gamma
  1. # σ or sigma the annualized volatility of the underlying stock
  2. eu_put = BinomialLRWithGreeks( 50, 52, r=0.05, T=2, N=2, sigma=0.3, is_put=True )
  3. results = eu_put.price()
  4.  
  5. print('European put values')
  6. print('Price: %s\nDelta: %s\nGamma: %s' % results)

 Nodes(stock prices): [array([84.90937139,              50. , 29.44315756]),
            array([113.61937767,  66.90626477,  39.3986339 ,  23.20040371]),
            array([152.03696329,  89.5289653 ,  52.72030863,  31.04504707,  18.28128424])]
the payoffs at maturity [ 0.          0.                 0.                     20.95495293   33.71871576]
At the time step 1 , payoffs: [ 0.          0.         10.39321697  26.26352636]
At the time step 0 , payoffs: [ 0.          5.15481754 17.75767426]
European put values
Price: 5.154817539524094
Delta: -0.3201529909525428
Gamma: 0.01678177314959919

VS 

     As we can see from the price() method and results, we managed to obtain additional information on Greeks from the modified binomial tree without any extra overhead in computational complexity. 

     Using the same example from the LR tree, we can compute the option values and Greeks for a European call and put option with 300 time steps:

  1. # σ or sigma the annualized volatility of the underlying stock
  2. eu_put = BinomialLRWithGreeks( 50, 52, r=0.05, T=2, N=300, sigma=0.3, 
  3.                                 is_put=True )
  4. results = eu_put.price()
  5.  
  6. print('European put values')
  7. print('Price: %s\nDelta: %s\nGamma: %s' % results)

 

  1. # σ or sigma the annualized volatility of the underlying stock
  2. eu_call = BinomialLRWithGreeks( 50, 52, r=0.05, T=2, N=300, sigma=0.3, 
  3.                                 is_put=False )
  4. results = eu_call.price()
  5.  
  6. print('European call values')
  7. print('Price: %s\nDelta: %s\nGamma: %s' % results)

     As we can see from the price() method and results, we managed to obtain additional information on Greeks from the modified binomial tree without any extra overhead in computational complexity.

Trinomial trees in option pricing

     In the binomial tree, each node leads to two other nodes in the next time step. Similarly, in a trinomial tree, each node leads to three other nodes in the next time step. Besides having up and down states, the middle node of the trinomial tree indicates no change in state. When extended over more than 2 time steps, the trinomial tree can be thought of as a recombining tree, where the middle nodes always retain the same values as the previous time step.
==>

     Let's consider the Boyle trinomial tree, where the tree is calibrated[ˈkælɪbreɪtɪd]校准使...标准化 so that the probability of up, down, and flat movements, u, d, and m with risk-neutral probabilities , , and are as follows:
     We can see that

  • recombines to m =1. With calibration, the no state movement m grows at a flat rate of 1 instead of at the risk-free rate.
  • The variable v is the annualized dividend yield,
  • and σ is the annualized volatility of the underlying stock.

     In general, with an increased number of nodes to process, a trinomial tree gives better accuracy than the binomial tree when fewer time steps are modeled, saving on computation speed and resources. The following diagram illustrates the stock price movements of a trinomial tree with 2 time steps:

A class for the trinomial tree option pricing model

     Let's create a TrinomialTreeOption class, inheriting from the BinomialTreeOption class. 

The methods for the TrinomialTreeOption are provided as follows:

  1. import math
  2. import numpy as np
  3.  
  4. class TrinomialTreeOption( BinomialTreeOption ):
  5.     def setup_parameters( self ):
  6.         # Calculate the jump sizes u, d and m
  7.         self.u = math.exp( self.sigma*math.sqrt( 2.*self.dt) )
  8.         self.m = 1
  9.         self.d = 1/self.u
  10.         
  11.         # Calculate the transition probabilities
  12.         # self.dt : self.T/float(self.N)
  13.         exponential_part = math.exp( self.sigma * math.sqrt(self.dt/2.) )
  14.         self.qu = ( ( math.exp( (self.r-self.div)*self.dt/2)
  15.                       - 1/exponential_part
  16.                     )/( exponential_part - 1/exponential_part )
  17.                   )**2
  18.         self.qd = ( ( exponential_part
  19.                       - math.exp( (self.r-self.div)*self.dt/2)
  20.                     )/( exponential_part - 1/exponential_part )
  21.                   )**2
  22.         
  23.         self.qm = 1 - self.qu - self.qd
  24.         
  25.     def init_stock_price_tree(self):
  26.         # Initialize a 2D tree at t=0
  27.         self.STs = [ np.array([self.S0]) ]
  28.         
  29.         for i in range( self.N ):
  30.             prev_nodes = self.STs[-1]
  31.             st = np.concatenate( ( prev_nodes*self.u,     # up state nodes
  32.                                    [prev_nodes[-1]*self.m,# middle nodes 
  33.                                     prev_nodes[-1]*self.d # down state nodes
  34.                                    ]
  35.                                  ) )
  36.             self.STs.append( st )    
  37.             
  38.     def traverse_tree( self, payoffs ):
  39.         """
  40.         Starting from the time the option expires, traverse
  41.         backwards and calculate discounted payoffs at each node
  42.         """
  43. #     # Returns the maximum payoff values between 
  44. #     # exercising the American option early 
  45. #     # and not exercising the option at all
  46. #     def check_early_exercise( self, payoffs, node ):
  47. #         if self.is_call:
  48. #             return np.maximum( payoffs, self.STs[node]-self.K )
  49. #         else:
  50. #             # payoffs: [1.41475309 9.46393007]
  51. #             # payoffs at maturity: [-8. 12.]
  52. #             # np.maximum ==>[1.41475309 12.]
  53. #             return np.maximum( payoffs, self.K-self.STs[node] )
  54.         
  55.         for i in reversed( range(self.N) ):
  56.             payoffs = ( payoffs[:-2] * self.qu +
  57.                         payoffs[1:-1] * self.qm + 
  58.                         payoffs[2:] * self.qd
  59.                       ) * self.df # df: The discount factor
  60.                       # * math.exp( -(self.r - self.div)*self.dt )   
  61.             if not self.is_european:
  62.                 payoffs = self.check_early_exercise( payoffs, i )
  63.                 
  64.             print("At the time step",i,", payoffs:",payoffs)
  65.                 
  66.         return payoffs

Using the same example of the binomial tree, we get the following result:

  1. eu_put = TrinomialTreeOption( 50, 52, r=0.05, T=2, N=2, sigma=0.3, is_put=True )
  2.  
  3. print( "European put option:", eu_put.price() )
Nodes: [array([50]), 
        array([76.42325802,   50.        , 32.71255459]), 
        array([116.81028732,  76.42325802, 50.        ,  32.71255459,  21.40222456])]
the payoffs at maturity [ 0.          0.          2.         19.28744541 30.59777544]
At the time step 1 , payoffs: [ 0.46353923  5.42131793 16.75137548]
At the time step 0 , payoffs: [6.57356527]
European put option: 6.573565269142496

  1. am_put = TrinomialTreeOption(50, 52, 
  2.                 r=0.05, T=2, N=2, sigma=0.3, is_put=True, is_am=True)
  3. print( "American put option:", am_put.price() )
Nodes: [array([50]), 
        array([76.42325802, 50.        , 32.71255459]), 
        array([116.81028732,76.42325802, 50.        ,  32.71255459, 21.40222456])]
the payoffs at maturity [ 0.          0.          2.         19.28744541 30.59777544]
At the time step 1 , payoffs: [ 0.46353923  5.42131793 19.28744541]
At the time step 0 , payoffs: [7.16134922]
American put option: 7.161349217272585

     By the trinomial tree model, we obtain prices of $6.57 and $7.16 for the European and American put options, respectively.

Lattices in option pricing

  • In binomial trees, each node recombines at every alternative可替代的 node.
  • In trinomial trees, each node recombines at every other node.

This property of recombining trees can also be represented as lattices[ˈlætɪs]格子,[计]点阵 to save memory without recomputing and storing recombined nodes.

Using a binomial lattice(CRR)

     We will create a binomial lattice from the binomial CRR tree(Cox-Ross-Rubinstein) since at every alternate up and down nodes, the prices recombine to the same probability of ud=1. (self.u : Expected value in the up state AND self.d : Expected value in the down state) In the following diagram,andrecombine with . The tree can now be represented as a single list:(
S*u*u = S_uu,
S_uu*d = S*u*u*d = S*u= S_u,
S_u*d = S*u*d = S = S_0
S*d = S_d
S_d * d = S_dd
 )
     For a N-step binomial tree, a list of size 2N +1 is required to contain the information on the underlying stock prices.

  • For European option pricing, the odd nodes of payoffs from the list represent the option value upon maturity. The tree traverses backward to obtain the option value.
  • For American option pricing, as the tree traverses backward, both ends of the list shrink列表的两端缩小, and the odd nodes represent the associated stock prices for any time step. Payoffs from the earlier exercise can then be taken into account.

A class for the CRR(Cox-Ross-Rubinstein) binomial lattice option pricing model 

     Let's convert the binomial tree pricing into a lattice by CRR. We can inherit from the BinomialCRROption class (which in turn inherits the BinomialTreeOption class) and create a new class named BinomialCRRLattice, as follows:

  1. import numpy as np
  2.  
  3. # import math
  4.  
  5. # # Price an option by the binomial CRR model
  6.  
  7. # class BinomialCRROption( BinomialTreeOption ):
  8. #     def setup_parameters(self):  # self.dt = self.T/float(self.N)
  9. #         # Expected stock price % in the up state  
  10. #         self.u = math.exp( self.sigma * math.sqrt(self.dt) )
  11. #         # Expected stock price % in the down state
  12. #         self.d = 1./self.u
  13. #
  14. #         # Calculate the transition probabilities
  15. #         # qu : The risk-neutral (up state) probability 
  16. #         self.qu = ( math.exp( (self.r-self.div)*self.dt )-self.d 
  17. #                   )/(self.u-self.d)
  18. #         # qd : The risk-neutral (down state) probability
  19. #         self.qd = 1-self.qu
  20.  
  21. class BinomialCRRLattice( BinomialCRROption ):
  22.     def setup_parameters( self ):
  23.         # inherit from the 'BinomialCRROption' class
  24.         # (which in turn inherits the 'BinomialTreeOption' class)
  25.         # python 2.x
  26.         super(BinomialCRRLattice, self).setup_parameters()
  27.         # python 3.x
  28.         # super().setup_parameters()
  29.         # N: number of time steps
  30.         self.M = 2*self.N + 1 # Number of nodes of Lattice
  31.     
  32.     def init_stock_price_tree( self ):
  33.         self.STs = np.zeros( self.M ) # initialize a lattice
  34.         # lattice[0]
  35.         self.STs[0] = self.S0 * self.u**self.N # e.g. if N=2, then S_uu=S*u*u
  36.         
  37.         for i in range(self.M)[1:]:            # u*d=1
  38.             self.STs[i] = self.STs[i-1]*self.d # if i=1, then S_uu*d=S*u*u*d=S*u
  39.                    # Current stock price <== S # if i=2, then S_u*d = S*u*d = S
  40.                                                # ...
  41.     def init_payoffs_tree( self ):
  42.         odd_nodes = self.STs[::2] # Take odd nodes only
  43.         print('the stock price at maturity',odd_nodes)
  44.         # calculate all payoffs, e.g. payoffs at S_uu, S_0, S_dd
  45.         if self.is_call: # call option
  46.             return np.maximum(0, odd_nodes-self.K)
  47.         else: # put option
  48.             return np.maximum(0, self.K-odd_nodes)
  49.         
  50.     def check_early_exercise( self, payoffs, node ):
  51.         # S_u, S_0, S_d
  52.         self.STs = self.STs[1:-1] # Shorten ends of the list
  53.         # S_uu, S_0, S_dd
  54.         odd_STs = self.STs[::2]   # Take odd nodes only
  55.         if self.is_call:
  56.             return np.maximum( payoffs, odd_STs - self.K )
  57.         else:
  58.             return np.maximum( payoffs, self.K - odd_STs )

     Using the same stock information from our binomial CRR model example, we can price a European and American put option using the binomial lattice pricing:

  1. eu_option = BinomialCRRLattice(
  2.                 50, 52, r=0.05, T=2, N=2, sigma=0.3, is_put=True)
  3.  
  4. print( 'European put option price:', eu_option.price() )

Nodes: [91.10594002 67.49294038 50. 37.04091103 27.4405818 ]
the stock price at maturity [91.10594002 50. 27.4405818 ] VS 52
the payoffs at maturity [ 0. 2. 24.5594182]
At the time step 1 , payoffs: [ 0.93269783 12.42301904]
At the time step 0 , payoffs: [6.24570845]
European put option price: 6.245708445206432
(save memory without recomputing and storing recombined nodes, e.g. Sud=Sdu=S0=50) VS The Cox–Ross–Rubinstein model for the code  def init_stock_price_tree(self):

  1. am_option = BinomialCRRLattice(50, 52,
  2.                 r=0.05, T=2, N=2, sigma=0.3, is_put=True, is_am=True)
  3.  
  4. print("American put option price:", am_option. price() )

Nodes: [91.10594002 67.49294038 50.         37.04091103 27.4405818 ]
the stock price at maturity [91.10594002 50.         27.4405818 ]
the payoffs at maturity [ 0.         2.        24.5594182]
At the time step 1 , payoffs: [ 0.93269783 12.42301904]
At the time step 1 , max payoffs: [ 0.93269783 14.95908897]
At the time step 0 , payoffs: [7.4284019]
At the time step 0 , max payoffs: [7.4284019]
American put option price: 7.428401902704828

VS The Cox–Ross–Rubinstein model for the code  def init_stock_price_tree(self):

     By using the CRR binomial tree lattice pricing model, we obtain prices of $6.2457 and $7.428 for the European and American put options, respectively.
https://www.maths.usyd.edu.au/u/UG/SM/MATH3075/r/Slides_7_Binomial_Market_Model.pdf

  1. #                           # delta_t=1/12
  2. print( "1+r*delta_t", 1+0.1*(1/12), 
  3.        "is close to", math.exp( 0.1*(1/12) )
  4.      )              
  5. #                    r=0.1  variance=0.1             2*sigma     
  6. print( "p~:", 0.5 + ( (0.1 - 0.1/2)*math.sqrt(1/12)/(2*math.sqrt(0.1) ) ) )
  7.  
  8. #                     sigma            T=4/12=1/3
  9. print( "u:",
  10.        math.exp( math.sqrt(0.1) * math.sqrt( (1/3)/float(4) ) )
  11.      )
  12.  
  13. print( "d:",
  14.        1/math.exp( math.sqrt(0.1) * math.sqrt( (1/3)/float(4) ) )
  15.      )

     For European option pricing, the odd nodes of payoffs from the list represent the option value upon maturity. The tree traverses backward to obtain the option value.

  1. eu_option2 = BinomialCRRLattice(50, 53,
  2.     r=0.1, T=4/12, N=4, sigma=math.sqrt(0.1), is_put=True, is_am=False)
  3.  
  4. print("European put option price is:", eu_option2.price())

(save memory without recomputing and storing recombined nodes) 

Nodes: [72.03638777 65.75161833 60.01515966 54.77917472 50.         45.63778138
 41.65614179 38.02187784 34.70468297]
the stock price at maturity [72.03638777 60.01515966 50.         41.65614179 34.70468297] vs 53
the payoffs at maturity [ 0.          0.          3.         11.34385821 18.29531703]
At the time step 3 , payoffs: [ 0.          1.4192295   6.92238713 14.53829067]
At the time step 2 , payoffs: [ 0.67140413  4.01086629 10.46784527]
At the time step 1 , payoffs: [2.24565342 7.03222631]
At the time step 0 , payoffs: [4.49143255]
European put option price is: 4.491432553637773

  1. am_option2 = BinomialCRRLattice(50, 53,
  2.     r=0.1, T=4/12, N=4, sigma=math.sqrt(0.1), is_put=True, is_am=True)
  3.  
  4. print("American put option price is:", am_option2.price())

Nodes: [72.03638777 65.75161833 60.01515966 54.77917472 50.         45.63778138
 41.65614179 38.02187784 34.70468297]
the stock price at maturity [72.03638777 60.01515966 50.         41.65614179 34.70468297]
the payoffs at maturity [ 0.          0.          3.         11.34385821 18.29531703]
At the time step 3 , payoffs: [ 0.          1.4192295   6.92238713 14.53829067]  #European
At the time step 3 , max payoffs: [ 0.          1.4192295   7.36221862 14.97812216] #American
##############def check_early_exercise( self, payoffs, node ):
self.STs = self.STs[1:-1] # Shorten ends of the list
Nodes: [72.03638777 65.75161833 60.01515966 54.77917472 50.         45.63778138
 41.65614179 38.02187784 34.70468297]
odd_STs = self.STs[::2]   # Take odd nodes only
Nodes: [ 65.75161833>53 60.01515966 54.77917472>53 50.         53-45.63778138=7.36221862
 41.65614179 53-38.02187784=14.97812216 ]
     For American option pricing, as the tree traverses backward, both ends of the list shrink列表的两端缩小, and the odd nodes represent the associated stock prices for any time step. Payoffs from the earlier exercise can then be taken into account.
##############
At the time step 2 , payoffs: [ 0.67140413  4.21894023 10.90402672]
At the time step 2 , max payoffs: [ 0.67140413  4.21894023 11.34385821]
At the time step 1 , payoffs: [2.34408831 7.55455975]
At the time step 1 , max payoffs: [2.34408831 7.55455975]
At the time step 0 , payoffs: [4.78958701]
At the time step 0 , max payoffs: [4.78958701]
American put option price is: 4.789587008231834
https://www.maths.usyd.edu.au/u/UG/SM/MATH3075/r/Slides_7_Binomial_Market_Model.pdf

Using the trinomial lattice

     The trinomial lattice works in very much the same way as the binomial lattice. Since each node recombines at every other node instead of alternate nodes, extracting odd nodes from the list is not necessary. Since the size of the list is the same as the one in the binomial lattice, there are no extra storage requirements in trinomial lattice pricing, as explained in the following diagram:

A class for the trinomial lattice option pricing model

     In Python, let's create a class named TrinomialLattice for the trinomial lattice implementation that inherits from the TrinomialTreeOption class(which in turn inherits the BinomialTreeOption class).

     Just as we did for the BinomialCRRLattice class, the setup_parameters, init_stock_price_tree , init_payoffs_tree , and check_early_exercise methods are overwritten, without having to take into account the payoffs at odd nodes:

  1. import numpy as np
  2.  
  3. # Price an option by the trinomial lattice
  4. class TrinomialLattice( TrinomialTreeOption ):
  5.     def setup_parameters( self ):
  6.         # inherit from the 'TrinomialTreeOption' class
  7.         # (which in turn inherits the 'BinomialTreeOption' class)
  8.         # python 2.x 
  9.         # super(TrinomialLattice, self).setup_parameters()
  10.         # python 3.x
  11.         super().setup_parameters()
  12.         self.M = 2*self.N + 1
  13.         
  14.     def init_stock_price_tree( self ):
  15.         self.STs = np.zeros( self.M ) # initialize a lattice
  16.         # lattice[0]
  17.         self.STs[0] = self.S0 * self.u**self.N # e.g. if N=2, then S_uu=S*u*u
  18.         
  19.         for i in range(self.M)[1:]:            # u*d=1
  20.             self.STs[i] = self.STs[i-1]*self.d # if i=1, then S_uu*d=S*u*u*d=S*u
  21.                    # Current stock price <== S # if i=2, then S_u*d = S*u*d = S
  22.                                                # ...
  23.     def init_payoffs_tree( self ):
  24.         # without having to take into account the payoffs at odd nodes ==>self.STs
  25.         
  26.         # calculate all payoffs, e.g. payoffs at S_uu, S_0, S_dd
  27.         if self.is_call: # call option
  28.             return np.maximum(0, self.STs-self.K)
  29.         else: # put option
  30.             return np.maximum(0, self.K-self.STs)
  31.         
  32.     def check_early_exercise( self, payoffs, node ):
  33.         # S_u, S_0, S_d
  34.         self.STs = self.STs[1:-1] # Shorten ends of the list
  35.         # without having to take into account the payoffs at odd nodes ==>self.STs
  36.         if self.is_call:
  37.             return np.maximum( payoffs, self.STs - self.K )
  38.         else:
  39.             return np.maximum( payoffs, self.K - self.STs )

     Using the same examples as before, we can price the European and American options using the trinomial lattice model:

  1. eu_option = TrinomialLattice(
  2.                 50, 52, r=0.05, T=2, N=2, sigma=0.3, is_put=True)
  3. print('European put option price:', eu_option.price() )

Nodes: [116.81028732  76.42325802  50.          32.71255459  21.40222456]
the payoffs at maturity [ 0.          0.          2.         19.28744541 30.59777544]
At the time step 1 , payoffs: [ 0.46353923  5.42131793 16.75137548]
At the time step 0 , payoffs: [6.57356527]
European put option price: 6.573565269142496
(save memory without recomputing and storing recombined nodes) VS

  1. am_option = TrinomialLattice(50, 52, 
  2.                 r=0.05, T=2, N=2, sigma=0.3, is_put=True, is_am=True)
  3. print('American put option price:', am_option.price() )

Nodes: [116.81028732  76.42325802  50.          32.71255459  21.40222456]
the payoffs at maturity [ 0.          0.          2.         19.28744541 30.59777544]
At the time step 1 , payoffs: [ 0.46353923  5.42131793 19.28744541]
At the time step 0 , payoffs: [7.16134922]
American put option price: 7.161349217272585 
(save memory without recomputing and storing recombined nodes) VS

     The output agrees with the results that were obtained from the trinomial tree option pricing model.

################################################################

Stochastics : https://blog.csdn.net/Linli522362242/article/details/96593981

Random Numbers

     It all starts with pseudo-random numbers, which build the basis for all simulation efforts; although quasi-random numbers准随机数 (e.g., based on Sobol sequences) have gained some popularity in finance, pseudo-random numbers伪随机数 still seem to be the benchmark.

     Throughout this chapter, to generate random numbers(For simplicity, we will speak of random numbers knowing that all numbers used will be pseudo-random), the functions provided by the numpy.random subpackage are used:

  1. import numpy as np
  2. # Imports the random number generation subpackage from NumPy.
  3. import numpy.random as npr
  4. import matplotlib.pyplot as plt
  5. %matplotlib inline

     For example, the rand() function returns random numbers from the open interval [0,1) in the shape provided as a parameter to the function. The return object is an ndarray object. Such numbers can be easily transformed to cover other intervals of the real line. For instance, if one wants to generate random numbers from the interval, one can transform the returned numbers from npr.rand() as in the next example—this also works in multiple dimensions due to NumPy broadcasting:

  1. npr.seed(100) # Fixes the seed value for reproducibility
  2. np.set_printoptions(precision=4) # fixes the number of digits for printouts
  3. # Uniformly distributed random numbers as one-dimensional ndarray object
  4. npr.rand(10)

 

npr.rand(5,5) # Uniformly distributed random numbers as 2-dimensional ndarray object.

 shape:(5,5)

  1. a = 5.  # Lower limit …
  2. b = 10. # … and upper limit …
  3. npr.rand(10) * (b - a) + a # for the transformation to another interval [5,10)

 

npr.rand(5, 5) * (b - a) + a # The same transformation [5,10) for 2 dimensions

 

  1. sample_size = 500
  2. rn1 = npr.rand(sample_size, 3)# Uniformly distributed random numbers in the given shape
  3. rn2 = npr.randint(0,10,sample_size)  # Random integers for a given interval [0,10)
  4. rn3 = npr.sample(size=sample_size)   # the “continuous uniform” distribution over the stated interval
  5. a = [0,25,50,75,100]
  6. rn4 = npr.choice(a, size=sample_size) # Randomly sampled values from a finite list object.
  7. rn4

 

  1. fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows = 2, ncols=2, figsize=(7,7))
  2.  
  3. ax1.hist(rn1, bins=25, stacked = True) #stacked by columns
  4. ax1.set_title('rand')
  5. ax1.set_ylabel('frequency')
  6. ax1.grid(True)
  7.  
  8. ax2.hist(rn2, bins=25)
  9. ax2.set_title('randint')
  10. ax2.grid(True)
  11.  
  12. ax3.hist(rn3, bins=25)
  13. ax3.set_title('sample')
  14. ax3.set_ylabel('frequency')
  15. ax3.grid(True)
  16.  
  17. ax4.hist(rn4, bins=25)
  18. ax4.set_title('choice')
  19. ax4.grid(True)

      Figure 12-1. Histograms of simple random numbers shows the results graphically for two continuous distributions(Uniformly distributed random numbers : rand, sample) and two discrete ones(randint(given interval [0,10)), choice(Randomly sampled values from a finite list object))

    Although there is much criticism[ˈkrɪtɪsɪzəm]批评 around the use of (standard) normal distributions in finance, they are an indispensable[ˌɪndɪˈspensəbl]不可或缺的 tool and still the most widely used type of distribution, in analytical as well as numerical applications. One reason is that many financial models directly rest in依赖于 one way or another on a normal distribution or a log-normal distribution. Another reason is that many financial models that do not rest directly on a (log-)normal assumption can be discretized离散化, and therewith对此 approximated for simulation purposes, by the use of the normal distribution. 

As an illustration, Figure 12-2 visualizes random draws from the following distributions:

  • Standard normal with mean of 0 and standard deviation of 1
  • Normal with mean of 100 and standard deviation of 20
  • Chi square with 0.5 degrees of freedom
  • Poisson with lambda of 1

     Figure 12-2 shows the results for the 3 continuous distributions and the discrete one (Poisson). The Poisson distribution is used, for example, to simulate the arrival of (rare) external events, like a jump in the price of an instrument or an exogenic shock外来冲击. Here is the code that generates it:

  1. sample_size = 500
  2. rn1 = npr.standard_normal(sample_size) # with mu=0, gamma=1
  3. rn2 = npr.normal(100,20, sample_size) # with mu=100 gamma=20  (standard deviation)
  4. rn3 = npr.chisquare(df=0.5, size=sample_size) # 0.5 degrees of freedom 
  5. rn4 = npr.poisson(lam=1.0, size=sample_size) #lambda =1 
  6.  
  7. fig, ( (ax1, ax2), (ax3,ax4) ) = plt.subplots(nrows=2, ncols=2, figsize=(7,7) )
  8.  
  9. ax1.hist(rn1, bins=25)
  10. ax1.set_title('standard normal')
  11. ax1.set_ylabel('frequency')
  12. ax1.grid(True)
  13.  
  14. ax2.hist(rn2, bins=25)
  15. ax2.set_title('normal(100, 20)')
  16. ax2.grid(True)
  17.  
  18. ax3.hist(rn3, bins=25)
  19. ax3.set_title('chi square')
  20. ax3.set_ylabel('frequency')
  21. ax3.grid(True)
  22.  
  23. ax4.hist(rn4, bins=25)
  24. ax4.set_title('Poisson')
  25. ax4.grid(True)

Figure 12-2. Histograms of random samples for different distributions 

Simulation

     In finance, 2 simulation tasks are of particular importance: simulation

  • of random variables and
  • of stochastic processes. 

     Monte Carlo simulation (MCS) is among the most important numerical techniques in finance, if not the most important and widely used. This mainly stems from the fact that it is the most flexible numerical method when it comes to the evaluation of mathematical expressions (e.g., integrals), and specifically the valuation of financial derivatives. The flexibility comes at the cost of a relatively high computational burden, though, since often hundreds of thousands or even millions of complex computations have to be carried out to come up with a single value estimate. 

Random Variables

     Consider, for example, the Black-Scholes-Merton setup for option pricing. In their setup, the level of a stock index at a future date given a level as of today stock index is given according to Equation 12-1.

Equation 12-1. Simulating future index level in Black-Scholes-Merton setup

The variables and parameters have the following meaning:

  •  : Index level at date 
  • : Constant riskless short rate
  • : Constant volatility (= standard deviation of returns) of
  • : Standard normally distributed random variable

     This financial model is parameterized and simulated as follows. The output of this simulation code is shown in Figure 12-3: 

  1. s0 = 100   # initial value
  2. r = 0.05   # constant short rate
  3. sigma=0.25 # constant volatility
  4. T=2.0      # in years # The horizon in year fractions
  5. Z=10000    # number of random draws
  6.  
  7. #the level of a stock index ST1 at a future date T given 
  8. # a level s0 as of today is given according to Equation
  9. # standard_normal(size) Samples from a standard normal distribution (mean=0, stdev=1)
  10.  
  11. ST1 = s0 * np.exp( (r-0.5 * sigma **2 )* T + sigma*np.sqrt(T) * npr.standard_normal(Z) )
  12.  
  13. plt.figure(figsize=(10,6))
  14. plt.hist(ST1, bins=50)
  15. plt.xlabel('index level')
  16. plt.ylabel('frequency')

 Figure 12-3. Statically simulated geometric Brownian motion (via npr.standard_normal())

     Figure 12-3 suggests that the distribution of the random variable as defined in Equation 12-1 is log-normal. One could therefore also try to use the numpy.random.lognormal() function to directly derive the values for the random variable. In that case, one has to provide the mean and the standard deviation to the function: 

  1. s0 = 100 # initial value
  2. r = 0.05 # constant short rate
  3. sigma=0.25 # constant volatility
  4. T=2.0 # in years
  5. Z=10000 # number of random draws
  6.  
  7. #the level of a stock index ST1 at a future date T given a level s0 as of today is given according to Equation
  8.  
  9. #ST1 = s0 * np.exp( (r-0.5 * sigma **2 )* T + sigma*np.sqrt(T) * npr.standard_normal(Z) )
  10.  
  11. ST2 = s0 * npr.lognormal( (r-0.5* sigma**2)*T, #mu==mean of lognormal Distribution
  12.                            sigma*np.sqrt(T),   #standard deviation of lognormal Distribution
  13.                            size=Z            
  14.                         )
  15. plt.hist(ST2, bins=50)
  16. plt.xlabel('index level')
  17. plt.ylabel('frequency')
  18. plt.grid(True)

 Figure 12-4. Statically simulated geometric Brownian motion (via npr.lognormal())

     By visual inspection, Figures 12-3 and 12-4 indeed look pretty similar. But let us verify this more rigorously by comparing statistical moments统计矩 of the resulting distributions.

     To compare the distributional characteristics of simulation results, the scipy.stats subpackage and the helper function print_statistics(), as defined here, prove useful

  1. import scipy.stats as scs
  2.  
  3. def print_statistics(a1, a2):
  4.     sta1 = scs.describe(a1)
  5.     sta2 = scs.describe(a2)
  6.     print('%14s %14s %14s' % ('statistic', 'data set 1', 'data set 2') )
  7.     print(45 * "-")
  8.     print('%14s %14.3f %14.3f' % ('size', sta1[0], sta2[0]))
  9.     print('%14s %14.3f %14.3f' % ('min', sta1[1][0], sta2[1][0]))
  10.     print('%14s %14.3f %14.3f' % ('max', sta1[1][1], sta2[1][1]))
  11.     print('%14s %14.3f %14.3f' % ('mean', sta1[2], sta2[2]))
  12.     print('%14s %14.3f %14.3f' % ('std', np.sqrt(sta1[3]), np.sqrt(sta2[3])))
  13.     print('%14s %14.3f %14.3f' % ('skew', sta1[4], sta2[4]))
  14.     print('%14s %14.3f %14.3f' % ('kurtosis', sta1[5], sta2[5]))
  15.  
  16. print_statistics(ST1, ST2)

 Obviously, the statistics of both simulation results are quite similar. The differences are mainly due to what is called the sampling error in simulation. Another error can also be introduced when discretely simulating continuous stochastic processes—namely the discretization error, which plays no role here due to the static nature of the simulation approach.

Stochastic Processes

     Roughly speaking, a stochastic process is a sequence of random variables. In that sense, we should expect something similar to a sequence of repeated simulations of a random variable when simulating a process随机变量的一系列 重复模拟. This is mainly true, apart from the fact that the draws are in general not independent but rather depend on the result(s) of the previous draw(s). In general, however, stochastic processes used in finance exhibit the Markov property, which mainly says that tomorrow’s value of the process only depends on today’s state of the process, and not any other more “historic” state or even the whole path history. The process then is also called memoryless.

###########################

Brownian motion(Wiener process维纳过程)

     以W(t）表示运动中一微粒从时刻t=0到时刻t>0的位移的横坐标（同样也可以讨论纵坐标），且设W(0)=0，根据爱因斯坦1905年提出的理论，微粒的这种运动是由于受到大量随机的相互独立的分子的碰撞的结果. 于是，粒子在时段（s,t]上的位移可以看作是许多微小位移的代数和. 则W(t)-W(s）服从正态分布。

The Wiener process is characterised by the following properties:

  • W has independent increments: for every t>0, the future increments, u>=0, are independent of the past values , s<=t,
    s <= t <= t+u
  • W has Gaussian increments: is normally distributed with mean 0 and variance u,  ~ 
  • W has continuous paths: is continuous in t.

     That the process has independent increments means that if then andare independent random variables, and the similar condition holds for n increments.

Basic properties of a one-dimensional Wiener process

     The unconditional probability density function follows a normal distribution with mean = 0 and variance = t, at a fixed time t: 

normal distribution:

The expectation is zero: 

The variance, using the computational formula, is t:  

These results follow immediately from the definition that increments have a normal distribution, centered at zero. Thus

stochastic differential equation (SDE) 

We assume that the stock price is driven by the stochastic differential equation (SDE)    (1)

     where is Brownian motion. We simulateover the time interval [0, T], which we assume to be is discretized as, where the time increments are equally spaced with width dt: Equally-spaced time increments is primarily used for notational convenience, because it allows us to writeas simply dt: All the results derived with equally-spaced increments are easily generalized to unequal spacing

Integrating from t to t + dt produces   (2)

     Equation (2) is the starting point for any discretization scheme. At time t, the value of is known, and we wish to obtain the next value  

Euler scheme

     The simplest way to discretize the process in Equation (2) is to use Euler discretization. This is equivalent to approximating the integrals using the leftpoint rule. Hence the first integral is approximated as the product of the integrand at time t, and the integration range dt

     We use the left-point rule since at time t the value is known. The right-hand rule would require that be known at time t. In an identical fashion, the second integral is approximated as

     since and are identical in distribution, where Z is a standard normal variable.(==). Hence, Euler discretization of (2)  is
   (3) 

Euler Scheme for the Black-Scholes Model 

    The Black-Scholes stock price dynamics under the risk neutral measure are  
                (4) OR 

 An application of Equation (3) produces Euler discretization for the BlackScholes model 

       (5)       : Constant riskless short rate

        ==> 
######

=
     A process S is said to follow a geometric Brownian motion with constant volatility σ and constant drift μ if it satisfies the stochastic differential equation OR , for a Brownian motion B. Applying Itô's lemma伊藤引理 with gives正常的全微分

， d_Y = d_St
==(6)

######

     Alternatively, we can generate log-stock prices, and exponentiate the result. By lemma ln follows the process 
   (6)

Euler discretization via Equation (3) produces 
so that  (7) OR 

where 

###########################

Geometric Brownian motion

     Consider now the Black-Scholes-Merton model in its dynamic form, as described by the stochastic differential equation (SDE, 随机微分方程) in Equation 12-2. Here, is a standard Brownian motion. The SDE is called a geometric Brownian motion. The values of are log-normally distributed and the (marginal) returns normally. 

Equation 12-2. Stochastic differential equation in Black-Scholes-Merton setup

     The SDE in Equation 12-2 can be discretized exactly by an Euler scheme. Such a scheme is presented in Equation 12-3, with being the fixed discretization interval and being a standard normally distributed random variable.

Equation 12-3. Simulating index levels dynamically in Black-Scholes-Merton setup 

  1. s0 = 100 # initial value
  2. T=2.0 # in years
  3. r = 0.05 # constant short rate
  4. sigma=0.25 # constant volatility
  5. Z=10000 # number of random draws
  6.  
  7. I =10000  #The number of paths to be simulated.
  8. M = 50    #The number of time intervals for the discretization.
  9. dt = T/M  #The length of the time interval in year fractions.
  10. S = np.zeros((M+1, I) )
  11. S[0] = s0 #first row #The initial values for the initial point in time t=0.
  12.  
  13.  
  14. #ST1 = s0 * np.exp( (r-0.5 * sigma **2 )* T + sigma*np.sqrt(T) * npr.standard_normal(Z) )
  15.  
  16. #ST2 = s0 * npr.lognormal( (r-0.5* sigma**2)*T, #mu==mean of lognormal Distribution
  17. #                           sigma*np.sqrt(T),   #standard deviation of lognormal Distribution
  18. #                           size=Z            
  19. #                        )
  20.  
  21. for t in range(1, M+1): #a stochastic process is a sequence of random variables
  22.     S[t] = S[t-1] * np.exp((r-0.5*sigma**2)*dt + sigma*np.sqrt(dt)*npr.standard_normal(Z) )
  23. # The simulation via semivectorized expression; 
  24. # the loop is over the points in time starting at t=1 and ending at t=T.
  25.  
  26. plt.figure( figsize=(10,6) )
  27. plt.hist(S[-1], bins=50) # the option price at maturity
  28. plt.xlabel('index level')
  29. plt.ylabel('frequency')

 Figure 12-5. Dynamically simulated geometric Brownian motion at maturity(S[-1])

     Following is a comparison of the statistics resulting from the dynamic simulation as well as from the static simulation(ST2).

print_statistics(S[-1], ST2)

 

Figure 12-6 shows the first 10 simulated paths:

  1. #Using the dynamic simulation approach allows us to visualize paths as displayed
  2. plt.figure(figsize=(10,6))
  3. plt.plot(S[:, :10], lw=1.5)
  4. plt.xlabel('time')
  5. plt.ylabel('index level')

 Figure 12-6. Dynamically simulated geometric Brownian motion paths

     Using the dynamic simulation approach not only allows us to visualize paths as displayed in Figure 12-6, but also to value估值

  • options with American/Bermudan exercise
  • or options whose payoff is path-dependent 收益与路径相关的期权. 

理解Self Financing 动态对冲

     Itô's lemma can be used to derive the Black–Scholes equation for an option. Suppose a stock price follows a geometric Brownian motion given by the stochastic differential equation dS = S(σdB + μ dt) <==. Then, if the value of an option at time t is f(t, St), Itô's lemma gives 

 

https://ethz.ch/content/dam/ethz/special-interest/mavt/dynamic-systems-n-control/idsc-dam/Lectures/Stochastic-Systems/SDE.pdf

 ...

 (26)==>https://www.zhihu.com/question/36828182?sort=created

     The term represents the change in (Option) value in time dt of the trading strategy consisting of holding an amount  (or)of the stock

 ################################################################ 

Finite differences in option pricing 

     Finite difference schemes有限差分方案 are very much similar to trinomial tree option pricing, where each node is dependent on three other nodes with an up movement, a down movement, and a flat movement. The motivation behind the finite differencing is the application of the Black-Scholes Partial Differential Equation (PDE 偏微分方程) framework (involving functions and their partial derivatives偏导数), where

  • price S(t) is a function of ƒ(S,t),
  • with r as the risk-free rate,
  • t as the time to maturity,
  • and σ as the volatility of the underlying security:
  • Black, Scholes and Merton showed that a riskless portfolio made up of an asset with value S and an option with value ƒ(S,t) satisfies the differential equation
  •  or 
    Equation 1: Black-Scholes-Merton PDE

    It is comprised of (partial) derivatives with repsect to time t and asset value S.

         The solution to the Black-Scholes-Merton PDE depends on several factors, including the expected form of ƒ(S,t) and boundary conditions imposed on the solution. Boundary conditions are specified to reflect the expected payoff of the option at expiry and for minimum and maximum values of S. 

     The finite difference technique tends to converge faster than lattices and approximates complex exotic options very well.

     To solve a PDE by finite differences working backward in time, a discrete-time grid of size M by N is set up to reflect asset prices over a course of time建立一个大小为 M×N 的离散时间网格以反映一段时间内的资产价格, so that S and t take on the following values at each point on the grid:
     The lattice is generated by dividing the time between today and expiry into M equal periods, and the underlying price into N equal levels. This leads to a grid with M+1 time points and N+1 price levels. The grid is typically chosen so that the current price of the underlying asset lies close to the middle of the N equal price levels of the grid.

     It follows that by grid notation,. is a suitably large asset price that cannot be reached by the maturity time, T. Thus dS and dt are intervals between each node in the grid, incremented by price and time, respectively.(Note that time is considered to run from 0 (today) to M*Δt=M*dt=T (at expiry). Alternative derivations may reverse this order and consider time counting down towards expiry.) 

The terminal condition at expiration time T for every value of S is

  • max(S − K, 0) for a call option with strike K and
  • max(K − S, 0) for a put option.

The grid traverses backward from the terminal conditions, complying with the PDE while adhering to the boundary conditions of the grid(遵守 PDE，同时遵守网格的边界条件), such as the payoff from an earlier exercise.

     The boundary conditions are defined values at the extreme ends of the nodes, where i=0 and i=N for every time at t. Values at the boundaries are used to calculate the values of all other lattice nodes iteratively using the PDE. 

     A visual representation of the grid is given in the following diagram. As i and j increase from the top-left corner of the grid, the price S tends toward (the maximum price possible) at the bottom-right corner of the grid:

Consider the function of one variable ƒ(x) shown in Figure 1.
     Then ƒ′(x), the derivative of the function at x, can be approximated in many ways. The most common are called the forward, backward and central approximation, all of which are drawn and indicated on Figure 1.

A number of ways to approximate the PDE are as follows: 

  • Forward difference:

    Consider the Taylor's series expansion for ƒ(x+h),
    Then re-arranging Equation 2 leads to,
         where Ο(h) symbolizes terms of the order of h( Ο(h) 表示 h 阶项 the truncation error). Assuming the step h is small then Ο(h) may be ignored and Equation 3 represents an approximation to ƒ′(x) at x.
  • Backward difference:
    Consider the Taylor's series expansion for ƒ(x-h),
    Then re-arranging Equation 4 leads to,
         where Ο(h) symbolizes terms of the order of h. Assuming the step h is small then Ο(h) may be ignored and Equation 5 represents an approximation to ƒ′(x) at x.
  • Central or symmetric difference:
         In addition to the forward and backwards approximation there is also a central approximation to ƒ′(x). This is formed by subtracting Equation 4 from Equation 2 and rearranging the result to obtain the equation,
         Note that the truncation error for the central approximation is Ο(h2), i.e. of the order of h2. This implies that it converges to the correct solution faster than either the forward or backward approximations.比前向或后向近似更快地收敛到正确的解决方案
  • The second derivative:
         In addition to the above first derivative approximations, to discretize the Black-Scholes-Merton PDE of Equation 1 also requires an approxmation for the second (partial) derivative term. This is formed by adding the two Talyor's series approximations of Equation 2and Equation 4and rearranging the result to obtain the equation,
         Note that as with the central approximation for the first derivative, the truncation error for this second derivative approximation is Ο(h2), i.e. of the order of h2. This has consequences for the speed of convergence of the finite difference algorithms.这对有限差分算法的收敛速度有影响。

     Once we have the boundary conditions set up, we can now apply an iterative approach using the explicit, implicit, or Crank-Nicolson method. 

The explicit method

The explicit method for approximatingƒ(S,t)==>is given by the following equation:
     Discretizing the Black-Scholes-Merton PDE : 
For the explicit method the Black-Scholes-Merton partial differential equation,
   is discretized using the following formulae:

where the indices i and j represent nodes on the pricing grid. 

Substituting these approximations into the PDE gives, 

     Here, we can see that the 1st difference is the backward difference with respect to t, the 2nd difference is the central difference with respect to S, and the 3rd difference is the second-order difference with respect to S. When we rearrange the terms, we get the following equation:

     By algebraic rearranging, we end up with the formula we need to implement to price the options, 
Where:
     
Then:

     The iterative approach of the explicit method can be visually represented by the following diagram:ƒ(S,t)

Put–call parity 

can be stated in a number of equivalent ways, most tersely as:
                 

     where 

  • C is the (current) value of a call, 
  • P is the (current) value of a put, 
  • D is the discount factor, 
  • F is the forward price of the asset,
  • and K is the strike price.
  • Note that the spot price is given by  (spot price现货价格 is present value, forward price is future value, discount factor relates these).
  • The left side corresponds to a portfolio of long a call and short a put,
    C - P : The minus sign indicates a short position.
  • while the right side corresponds to a forward contract.
  • The assets C and P on the left side are given in current values, while the assets F and K are given in future values (forward price of asset, and strike price paid at expiry), which the discount factor D converts to present values.

     Using spot price S instead of forward price F yields: 
 <== <==

Rearranging the terms yields a different interpretation: 

     In this case the left-hand side is a fiduciary call, which is long a call and enough cash (or bonds) to pay the strike price if the call is exercised,

     while the right-hand side is a protective put, which is long a put and the asset, so the asset can be sold for the strike price if the spot is below strike at expiry.

     Both sides have payoff max(S(T), K) at expiry (i.e., at least the strike price, or the value of the asset if more), which gives another way of proving or interpreting put–call parity. 

In more detail, this original equation can be stated as:
     

  • C(t) is the value of the call at time t,
  • P(t) is the value of the put of the same expiration date,
  • S(t) is the spot price of the underlying asset,
  • K is the strike price, and
  • B(t,T) is the present value of a zero-coupon bond that matures to $1 at time T. This is the present value factor for K.

     Note that the right-hand side of the equation is also the price of buying a forward contract on the stock with delivery price K. Thus one way to read the equation is that a portfolio that is long a call and short a put is the same as being long a forward. In particular, if the underlying is not tradeable but there exists forwards on it, we can replace the right-hand-side expression by the price of a forward. 

     If the bond interest rate, r, is assumed to be constant then 

     Note: r refers to the force of interest, which is approximately equal to the effective annual rate for small interest rates. However, one should take care with the approximation, especially with larger rates and larger time periods. To find r exactly, use r=ln(1+i), where i is the effective annual interest rate. 

     When valuing European options written on stocks with known dividends that will be paid out during the life of the option, the formula becomes:

where D(t) represents the total value of the dividends from one stock share to be paid out over the remaining life of the options, discounted to present value. We can rewrite the equation as:

and note that the right-hand side is the price of a forward contract on the stock with delivery price K, as before. 

####################################

     A stock price is currently $50. Over each of the next two(n_t=2) three-month( t=3/12, T=6/12) years) periods it is expected to go up by 6%(u=1+0.06=1.06) or down by 5%(d=1-0.05=0.95). The risk-free interest rate is 5% per annum with continuous compounding. What is the value of a six-month European call option with a strike price of $51?


Figure S12.1

     A tree describing the behavior of the stock price is shown in Figure S12.1. The risk-neutral probability of an up move, p(), is given by

     There is a payoff from the option of 56.18-51=5.18 ( <== max(S − K, 0) )  for the highest final node (which corresponds to two up moves) zero in all other cases. The value of the option is therefore

     For the situation considered in Problem 12.12, what is the value of a six-month European put option with a strike price of $51? Verify that the European call and European put prices satisfy put–call parity. If the put option were American, would it ever be optimal to exercise it early at any of the nodes on the tree?
Figure S12.2
     The tree for valuing the put option is shown in Figure S12.2. We get
a payoff of 51-50.35=0.65 (<==max(K − S, 0) ) if the middle final node is reached
and a payoff of 51-45.125=5.875 (<==max(K − S, 0) ) if the lowest final node is reached. The value of the option is therefore

 

The value of the put plus the stock price is from Problem 12.12

=

 The value of the call plus the present value of the strike price is

 ==> the discount factor D===>

 This verifies that put–call parity holds

####################################

Writing the finite difference base class 

     Since we will be writing the explicit, implicit, and Crank-Nicolson methods of finite differences in Python, let's write a base class that inherits the common properties and functions of all three methods. 

     We will create a class called FiniteDifferences that accepts and assigns all the required parameters in the __init__ constructor method. The price() method is the entry point for invoking the specific finite difference scheme implementation, and will invoke these methods in the following order: setup_boundary_conditions(), setup_coefficients(), traverse_grid(), and interpolate(). These methods are explained as follows:

  • setup_boundary_conditions : Sets up the boundary conditions of the grid structure as a NumPy two-dimensional array

    Dirichlet Condition
         A Dirichlet boundary condition is a boundary condition which assigns a value to f at boundary. For example, for a call option at S_min=0 , the option value could be set to zero, since the option is worthless. However, if underlying price reaches to S_max at time t, the value of option is the discounted value of its expected payoff.

          For a vanilla call option: ƒ(S=0,t) = 0
                                                    The boundary equation states that on the left-boundary, for all times up to expiry, the value of the call will equal 0.
                                                  ƒ(S_max,t) = S_max -
                                                  ƒ(S_max,t) = S_max - K * e^(-r * dt*(N - time_index) )
                                                  ,The equation for the right-boundary states that, for all times up to expiry, the value of the call will equal the pay-off function, albeit slightly adjusted to take into account discounting due to the risk-free rate, r. This is a consequence of Put-Call Parity.andand as S→∞ (S >>K), for all t,


          For a vanilla put option: ƒ(S=0,t) = , S_min=0 ==> C(0,t)=0
                                                  ƒ(S=0,t) = K * e^(-r * dt*(N - time_index) )
                                                  ƒ(S_max,t) = 0（当资产价格S-->S_max非常大时，put option几乎是一文不值，因为我们几乎可以肯定的是Strike price(K)低于现价）
           
  • setup_coefficients : Sets up the necessary coefficients that are used for traversing the grid structure
  • traverse_grid : Iterates the grid structure backward in time, storing the calculated values toward the first column of the grid
  • interpolate : Using the final calculated values on the first column of the grid, this method will interpolate these values to find the option price that closely infers the initial stock price, S0

     All of these methods are abstract methods that can be implemented by the derived classes. An exception type of NotImplementedError will be thrown should we forget to implement these methods.

The base class with the mandatory methods should look like this:

https://www.semanticscholar.org/paper/Excel-implementation-of-finite-difference-methods-Kyng-Purcal/cd18bc626774407cc692e562f1f4cbe7143942ca/figure/5

     Abstract base classes (ABCs) provide a way to define interfaces for a class. The @abstractmethod() decorator declares abstract methods that child classes should implement. Unlike Java's abstract methods, these methods may have an implementation and may be called via the super() mechanism from the class that overrides it.

     In addition to these methods, we would need to define dS and dt , the change in S per unit time, and the change in T per iteration, respectively. 

     Finally, add the price() method as the entry point that shows the steps in calling our discussed abstract methods:

  1. from abc import ABCMeta, abstractmethod
  2. import numpy as np
  3.  
  4. """ 
  5. Base class for sharing attributes and functions of FD 
  6. """
  7.  
  8. class FiniteDifferences( object ):
  9.     def __init__( self, S0, K, r=0.05, T=1,
  10.                   sigma=0, Smax=1, M=1, N=1, is_put=False
  11.                 ):
  12.         self.S0 = S0
  13.         self.K = K
  14.         self.r = r
  15.         self.T = T
  16.         self.sigma = sigma
  17.         self.Smax = Smax       
  18.         # N: number of time steps
  19.         self.M, self.N = M, N # stock price levels X time points
  20.         self.is_call = not is_put
  21.         
  22.         self.i_values = np.arange( self.M ) # Smin~Smax: price level indices
  23.         self.j_values = np.arange( self.N ) # 0~T: time point indices
  24.         
  25.         self.boundary_conds = np.linspace( 0, Smax, self.M+1 )
  26.         self.grid = np.zeros( shape=(self.M+1, self.N+1) )
  27.         """
  28.         https://numpy.org/doc/stable/reference/generated/numpy.linspace.html
  29.         numpy.linspace( start, stop, num=50, endpoint=True,
  30.                         retstep=False, dtype=None, axis=0 )
  31.         np.linspace(2.0, 3.0, num=5)
  32.         array([2.  , 2.25, 2.5 , 2.75, 3.  ])
  33.         """
  34.     @property
  35.     def dS( self ):
  36.         return self.Smax/float(self.M)
  37.     
  38.     @property
  39.     def dt( self ):
  40.         return self.T/float(self.N)
  41.     
  42.     @abstractmethod
  43.     def setup_boundary_conditions( self ):
  44.         raise NotImplementedError('Implementation required: boundary_conditions!')
  45.         
  46.     @abstractmethod
  47.     def setup_coefficients( self ):
  48.         raise NotImplementedError('Implementation required: coefficients!')
  49.     
  50.     @abstractmethod
  51.     def traverse_grid( self ):
  52.         # Iterate the grid backwards in time
  53.         raise NotImplementedError('Implementation required: traverse_grid!')
  54.     
  55.     @abstractmethod
  56.     def interpolate( self ):
  57.         """
  58.             Use piecewise linear interpolation on the initial
  59.             grid column to get the closest price at S0.
  60.             
  61.         numpy.interp(x, xp, fp, left=None, right=None, period=None)
  62.         left: Value to return for x < xp[0], default is fp[0]
  63.         right: Value to return for x > xp[-1], default is fp[-1]
  64.         One-dimensional linear interpolation for monotonically 
  65.         increasing sample points.
  66.         Returns the one-dimensional piecewise linear interpolant to
  67.         a function with given discrete data points (xp, fp), 
  68.         evaluated at x.
  69.         xp = [1, 2, 3]
  70.         fp = [3, 2, 0]
  71.         np.interp(2.5, xp, fp) ==> 1.0
  72.         
  73.         self.boundary_conds = np.linspace( 0, Smax, self.M+1 )
  74.         # grid : M+1 stock price levels X N+1 time points
  75.         """
  76.         # self.boundary_conds = np.linspace( 0, Smax, self.M+1 )
  77.         # self.grid = np.zeros( shape=(self.M+1, self.N+1) )
  78.         return np.interp( self.S0,
  79.                           self.boundary_conds, self.grid[:,0]
  80.                         )# lookup grid to return option present value
  81.     def price( self ):
  82.         self.setup_boundary_conditions()
  83.         self.setup_coefficients()
  84.         self.traverse_grid()
  85.         return self.interpolate()

pricing European options using the explicit method of finite differences

     The Python implementation of finite differences by using the explicit method is given in the following FDExplicitEu class, which inherits from the FiniteDifferences class and overrides the required implementation methods (with @abstractmethod() decorator):

  1. import numpy as np
  2.  
  3. # Explicit method of Finite Differences
  4.  
  5. class FDExplicitEu( FiniteDifferences ):
  6.     def setup_boundary_conditions(self):
  7.         print('The boundary(Stock price) at expiry:', self.boundary_conds)
  8.         print('strike price K:',self.K)
  9.         
  10.         if self.is_call:
  11.             # The terminal condition at expiration time T for every value of S :
  12.             # max(S − K, 0) for a call option
  13.             # max(K − S, 0) for a put option
  14.             # self.boundary_conds = np.linspace( 0, Smax, self.M+1 )
  15.             self.grid[:,-1] = np.maximum( 0, self.boundary_conds-self.K )
  16.             print("Payoff=max(boundary_conds-K,0):\n", self.grid )
  17.             
  18.             # N: number of time steps
  19.             # self.j_values = np.arange( self.N ) # 0~T: time point indices
  20.             # self.dt = self.T/float(self.N)
  21.             # grid: (M+1) stock price levels X (N+1) time points
  22.             # 1st time step: payoff  * 
  23.             self.grid[-1, :-1] = self.Smax-self.K*\
  24.                 np.exp( -self.r*self.dt*(self.N-self.j_values) )
  25.             print('Dirichlet boundary condition for call option',
  26.                   self.grid 
  27.                  )
  28.             
  29.         else:
  30.             self.grid[:,-1] = np.maximum( 0, self.K-self.boundary_conds )
  31.             print("Payoff=max(K-boundary_conds,0):\n", self.grid )
  32.             
  33.             self.grid[0, :-1] = self.K*\
  34.                 np.exp( -self.r*self.dt*(self.N-self.j_values) )
  35.             print('Dirichlet boundary condition for put option',
  36.                    self.grid 
  37.                  )       
  38.             
  39.     def setup_coefficients(self):
  40.         # self.i_values = np.arange( self.M ) # Smin~Smax: price level indices                        
  41.         self.a = 0.5*self.dt*( (self.sigma**2)*(self.i_values**2)-\
  42.                                 self.r*self.i_values
  43.                              )
  44.                                  
  45.         self.b = 1 - self.dt*( (self.sigma**2)*(self.i_values**2)+\
  46.                                 self.r
  47.                              )
  48.         self.c = 0.5*self.dt*( (self.sigma**2)*(self.i_values**2)+\
  49.                                 self.r*self.i_values
  50.                              )
  51.     def traverse_grid(self):
  52.          for j in reversed(self.j_values):
  53.              for i in range(self.M)[2:]:# i>=2 since S_min=S_0 = 0
  54.                  self.grid[i,j] = \
  55.                      self.a[i]*self.grid[i-1,j+1] +\
  56.                      self.b[i]*self.grid[i,j+1] + \
  57.                      self.c[i]*self.grid[i+1,j+1]
  58.          print( 'Option price grid', self.grid)

    On completion of traversing the grid structure, the 1st column contains the present value of the initial asset prices at t=0. The interp function of NumPy is used to perform a linear interpolation to approximate the option value.

     Besides using linear interpolation as the most common choice for the interpolation method, the other methods such as the spline or cubic may be used to approximate the option value.

  1. option = FDExplicitEu(30, 30, r=0.1, T=0.75, 
  2.     sigma=0.4, Smax=60, M=6, N=3, is_put=True)
  3. print(option.price())

The boundary(Stock price) at expiry: [ 0. 10. 20. 30. 40. 50. 60.]
          self.boundary_conds = np.linspace( 0, Smax, self.M+1 ), t=0
strike price K: 30
          self.grid[:,-1] = np.maximum( 0, self.K-self.boundary_conds )
Payoff=max(K-boundary_conds,0):
 [[ 0.  0.  0. 30.]
  [ 0.  0.  0. 20.]
  [ 0.  0.  0. 10.]
  [ 0.  0.  0.  0.]
  [ 0.  0.  0.  0.]
  [ 0.  0.  0.  0.]
  [ 0.  0.  0.  0.]
]
Dirichlet boundary condition for put option
          self.grid[0, :-1] = self.K*np.exp( -self.r*self.dt*(self.N-self.j_values) )
 [[27.8323 28.5369 29.2593 30.    ]
  [ 0.          0.           0.           20.    ]
  [ 0.          0.           0.           10.    ]
  [ 0.          0.           0.             0.    ]
  [ 0.          0.           0.             0.    ]
  [ 0.          0.           0.             0.    ]
  [ 0.          0.           0.             0.    ]
]
Option price grid
[[27.8323 28.5369 29.2593 30.    ]
 [ 0.            0.          0.       20.    ]
 [ 6.4964    7.6884  9.25   10.    ]
 [ 2.5289    2.1945  1.425   0.    ]
 [ 0.7214    0.3848  0.         0.    ]
 [ 0.1683    0.          0.         0.    ]
 [ 0.            0.          0.         0.    ]
]
option. price() ==>put option price: 2.5288940625

  1. option = FDExplicitEu(30, 30, r=0.1, T=0.75, 
  2.                       sigma=0.4, Smax=60, M=6, N=3, is_put=False)
  3. print("call option price:",option.price())

The boundary(Stock price) at expiry: [ 0. 10. 20. 30. 40. 50. 60.]
                         self.boundary_conds = np.linspace( 0, Smax, self.M+1 ), t=0
strike price K: 30
                         self.grid[:,-1] = np.maximum( 0, self.boundary_conds-self.K )
Payoff=max(boundary_conds-K,0):
 [[ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0.  0.  0. 10.]
 [ 0.  0.  0. 20.]
 [ 0.  0.  0. 30.]]
Dirichlet boundary condition for call option
[[ 0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.    ]
 [ 0.      0.      0.     10.    ]
 [ 0.      0.      0.     20.    ]
 [32.1677 31.4631 30.7407 30.    ]]                
                    self.grid[-1, :-1] = self.Smax-self.K* np.exp( -self.r*self.dt*(self.N-self.j_values) )
Option price grid
[[ 0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.    ]
 [ 0.5721  0.2284  0.      0.    ]
 [ 4.874   3.6758  2.175   0.    ]
 [12.9137 11.866  10.75   10.    ]
 [22.3525 21.476  20.75   20.    ]
 [32.1677 31.4631 30.7407 30.    ]]
call option price: 4.873984687 

     Consider the example of a European put option. The underlying stock price is $50 with a volatility of 40 percent(annualized volatility of the underlying stock : sigma=0.4). The strike price of the put option is $50 with an expiration time of 5 months. The risk-free rate is 10 percent. 

  1. option = FDExplicitEu(50, 50, r=0.1, T=5./12., 
  2.     sigma=0.4, Smax=100, M=100, N=1000, is_put=True)
  3. print(option.price())

 

     What happens when other values of M(stock price levels) and N(number of time steps) are chosen improperly?

  1. option = FDExplicitEu(50, 50, r=0.1, T=5./12., 
  2.                      sigma=0.4, Smax=100, M=80, N=100, is_put=True)
  3. print(option.price())

 

     when other values of M and N are chosen improperly, it appears that the explicit method of the finite difference scheme suffers from instability problems.

The implicit method

     The instability problem of the explicit method can be overcome using the forward difference with respect to time. The implicit method for approximatingƒ(S,t)==>  is given by the following equation:
For the implicit method the Black-Scholes-Merton partial differential equation,​​​​​​​
is discretized using the following formulae

  • use a forward approximation for ∂ƒ/∂t  (Compare this with the explicit method where the backward approximation is used.)
  • use a central approximation for ∂ƒ/∂S : 
  • use a standard approximation for : 

where the indices i and j represent nodes on the pricing grid. 

Substituting these approximations into the PDE gives, 

     Here, it can be seen that the only difference between the implicit and explicit approximating scheme lies in the first difference, where the forward difference with respect to t is used in the implicit scheme. When we rearrange the terms, we get the following expression: ​​​​​​​​​​​​​​​​​​​​​     
Where:

Here: 

     The iterative approach of the implicit scheme can be visually represented with the following diagram: ƒ(S,t)

     From the preceding diagram, we can note that values of j+1 need to be computed before they can be used in the next iterative step, as the grid traverses backward. In the implicit scheme, the grid can be thought of as representing a system of linear equations at each iteration, as follows: 
By rearranging the terms, we get the following equation:

     The linear system of equations can be represented in the form of Ax = B, where we want to solve values of x in each iteration. Since the matrix A is tri-diagonal三对角矩阵, we can use the LU factorization, where A=LU, for faster computation. Remember that we solved the linear system of equations using LU decomposition(or also known as lower-upper factorization) in https://blog.csdn.net/Linli522362242/article/details/125546725, The Importance of Linearity in Finance. 

     The LU decomposition, or also known as lower-upper factorization, is one of the methods that solve square systems of linear equations. As its name implies, the LU factorization decomposes the A matrix into a product of two matrices: a lower triangular matrix, L, and an upper triangular matrix, U. The decomposition can be represented as follows: 
We will look at a simple example and write out the row operations.

L : row operation and constraint : column_index= row_index
                               row_index+1=3               row_index+1=2https://blog.csdn.net/Linli522362242/article/details/125546725


==>(==>)==> vector y

==>(==>)==>vector x 

pricing European options using the implicit method of finite differences

     The Python implementation of the implicit scheme is given in the following FDImplicitEu class. We can inherit the implementation of the explicit method from the FDExplicitEu class we discussed earlier and override the necessary methods of interest, namely, the setup_coefficients and traverse_grid methods: 

  1. import numpy as np
  2. import scipy.linalg as linalg
  3.  
  4. # Explicit method of Finite Differences
  5.  
  6. class FDImplicitEu( FDExplicitEu ):
  7.     def setup_coefficients( self ):
  8.         # self.i_values = np.arange( self.M ) # Smin~Smax: price level indices
  9.         self.a = 0.5*self.dt*( self.r*self.i_values -\
  10.                                (self.sigma**2)*(self.i_values**2)
  11.                              )
  12.                                  
  13.         self.b = 1 + self.dt*( self.r +\
  14.                                (self.sigma**2)*(self.i_values**2)
  15.                              )
  16.         self.c =-0.5*self.dt*( self.r*self.i_values +\
  17.                                (self.sigma**2)*(self.i_values**2)
  18.                              )
  19.         # since S_min=S_0 = grid[0] : Boundary conditions
  20.         # grid: (M+1) stock price levels X (N+1) time points
  21.         self.coeffs = np.diag( self.a[2:self.M], 
  22.                                -1 #  k<0 for diagonals below the main diagonal.
  23.                              )+\
  24.                       np.diag( self.b[1:self.M],
  25.                                0  # default k=0 for diagonals
  26.                              )+\
  27.                       np.diag( self.c[1:self.M-1],
  28.                                1  # k>0 for diagonals above the main diagonal
  29.                              )
  30.     def traverse_grid( self ):
  31.         """ Solve using linear systems of equations """
  32.         P, L, U = linalg.lu( self.coeffs )
  33.         
  34.         aux = np.zeros( self.M-1 ) # 1~M-1
  35.         
  36.         for j in reversed( range(self.N) ): # 0~t~N-1 ==> N-1~t~0
  37.             aux[0] = np.dot( -self.a[1], self.grid[0,j] )
  38.             aux[-1] = np.dot( -self.c[-1], self.grid[-1, j] )
  39.             
  40.             # L@y = B ==>y
  41.             y = linalg.solve( L, 
  42.                               self.grid[1:self.M, j+1] + aux 
  43.                             )
  44.             # U@x = y ==>x
  45.             x = linalg.solve( U,
  46.                               y
  47.                             )
  48.             self.grid[1:self.M, j] = x
  49.         print( 'Option price grid', self.grid)

 

  1. option = FDImplicitEu(30, 30, r=0.1, T=0.75, 
  2.                       sigma=0.4, Smax=60, M=6, N=3, is_put=True)
  3. print("put option price:",option.price())

The boundary(Stock price) at expiry: [ 0. 10. 20. 30. 40. 50. 60.]
          self.boundary_conds = np.linspace( 0, Smax, self.M+1 ), t=0
strike price K: 30
          self.grid[:,-1] = np.maximum( 0, self.K-self.boundary_conds )
Payoff=max(K-boundary_conds,0):
 [[ 0.  0.  0. 30.]
  [ 0.  0.  0. 20.]
  [ 0.  0.  0. 10.]
  [ 0.  0.  0.  0.]
  [ 0.  0.  0.  0.]
  [ 0.  0.  0.  0.]
  [ 0.  0.  0.  0.]
]
Dirichlet boundary condition for put option
          self.grid[0, :-1] = self.K*np.exp( -self.r*self.dt*(self.N-self.j_values) )
 [[27.8323 28.5369 29.2593 30.    ]
  [ 0.          0.           0.           20.    ]
  [ 0.          0.           0.           10.    ]
  [ 0.          0.           0.             0.    ]
  [ 0.          0.           0.             0.    ]
  [ 0.          0.           0.             0.    ]
  [ 0.          0.           0.             0.    ]
]
Option price grid
[[27.8323 28.5369 29.2593 30.    ]
 [17.8959 18.5711 19.2729 20.    ]
 [ 8.5912  8.9637  9.4216 10.    ]
 [ 2.2156  1.7045  0.996   0.    ]
 [ 0.6429  0.4016  0.1697  0.    ]
 [ 0.1907  0.1049  0.0367  0.    ]
 [ 0.      0.      0.      0.    ]]
put option price: 2.2155510303909964

  1. print('coefficients:')
  2. print('a:',option.a)
  3. print('b:',option.b)
  4. print('c:',option.c)

  1. option = FDImplicitEu(30, 30, r=0.1, T=0.75, 
  2.                       sigma=0.4, Smax=60, M=6, N=3, is_put=False)
  3. print("call option price:",option.price())

The boundary(Stock price) at expiry: [ 0. 10. 20. 30. 40. 50. 60.]
                         self.boundary_conds = np.linspace( 0, Smax, self.M+1 ), t=0
strike price K: 30
                         self.grid[:,-1] = np.maximum( 0, self.boundary_conds-self.K )
Payoff=max(boundary_conds-K,0):
 [[ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0.  0.  0. 10.]
 [ 0.  0.  0. 20.]
 [ 0.  0.  0. 30.]]
Dirichlet boundary condition for call option
[[ 0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.    ]
 [ 0.      0.      0.     10.    ]
 [ 0.      0.      0.     20.    ]
 [32.1677 31.4631 30.7407 30.    ]]                
                    self.grid[-1, :-1] = self.Smax-self.K* np.exp( -self.r*self.dt*(self.N-self.j_values) )

Option price grid
[[0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00]
 [3.8233e-02 1.6884e-02 4.6786e-03 0.0000e+00]
 [7.3338e-01 4.0931e-01 1.5331e-01 0.0000e+00]
 [4.3584e+00 3.1504e+00 1.7278e+00 0.0000e+00]
 [1.2789e+01 1.1849e+01 1.0902e+01 1.0000e+01]
 [2.2344e+01 2.1557e+01 2.0771e+01 2.0000e+01]
 [3.2168e+01 3.1463e+01 3.0741e+01 3.0000e+01]]
call option price: 4.3584400028969705

     Using the same example as the explicit scheme, we can price the European put options using the implicit scheme: 

  1. option = FDImplicitEu(50, 50, r=0.1, T=5./12., 
  2.     sigma=0.4, Smax=100, M=100, N=1000, is_put=True)
  3. print(option.price())

  1. option = FDImplicitEu(50, 50, r=0.1, T=5./12., 
  2.     sigma=0.4, Smax=100, M=80, N=100, is_put=True)
  3. print(option.price())

 

     Given the current parameters and input data, we can see that there are no stability issues with the implicit scheme. 

  1. option = FDImplicitEu(140, 150, r=0.04, T=6/12, 
  2.                       sigma=0.25, Smax=200, M=200, N=200, is_put=False)
  3. print("call option price:",option.price())

  1. option = FDImplicitEu(140, 150, r=0.04, T=6/12, 
  2.                       sigma=0.25, Smax=200, M=200, N=200, is_put=True)
  3. print("put option price:",option.price())

 

The Crank-Nicolson method

     Another way of avoiding the instability issue, as seen in the explicit method, is to use the Crank-Nicolson method. The Crank-Nicolson method converges much more quickly using a combination of
the explicitand implicitmethods, taking the average of both. This leads us to the following equation:
and
==>j向左偏移1个单位==>

Crank-Nicolson method: 

 
\large +
                                                          
                                                         

This equation can also be rewritten as follows: 

\LARGE -\alpha_{i}f_{i-1,j-1} + (1-\beta_{i})f_{i,j-1} - \gamma_{i}f_{i+1,j-1}=\alpha_{i}f_{i-1,j} + (1+\beta_{i})f_{i,j} + \gamma_{i}f_{i+1,j}

Where: OR
     The iterative approach of the implicit scheme can be visually represented with the following diagram:vs

 We can treat the equations as a system of linear equations in a matrix form: https://www.goddardconsulting.ca/option-pricing-finite-diff-crank-nicolson.htmlhttps://quintus-zhang.github.io/post/on_pricing_options_with_finite_difference_methods/==> f_0=0 then j+1 ==> M_1 @ f_j = M_2 @ f_j+1

Where: 

We can solve for the matrix M on every iterative procedure. 

pricing European options using the Crank-Nicolson method of finite differences

     The Python implementation of the Crank-Nicolson method is given in the following FDCnEu class, which inherits from the FDExplicitEu class and overrides only the setup_coefficients and traverse_grid methods:​​​​​​​

  1. import numpy as np
  2. import scipy.linalg as linalg
  3.  
  4. # Crank-Nicolson method of Finite Differences
  5. class FDCnEu( FDExplicitEu ):
  6.     def setup_coefficients(self):
  7.         self.alpha = 0.25 * self.dt * ( (self.sigma**2)*(self.i_values**2) -\
  8.                                         self.r * self.i_values
  9.                                       )
  10.         self.beta = -0.5 * self.dt * ( (self.sigma**2)*(self.i_values**2) +\
  11.                                         self.r
  12.                                     )
  13.         self.gamma = 0.25 * self.dt * ( (self.sigma**2)*(self.i_values**2) +\
  14.                                         self.r * self.i_values
  15.                                       )
  16.         # since S_min=S_0 = grid[0] : Boundary conditions
  17.         # grid: (M+1) stock price levels X (N+1) time points
  18.         self.M1 = -np.diag( self.alpha[2:self.M], 
  19.                             -1 #  k<0 for diagonals below the main diagonal.
  20.                           )+\
  21.                   np.diag( 1-self.beta[1:self.M],
  22.                           0 # default k=0 for diagonals
  23.                          )-\
  24.                   np.diag( self.gamma[1:self.M-1],
  25.                           1 # k>0 for diagonals above the main diagonal
  26.                          )
  27.         self.M2 = np.diag( self.alpha[2:self.M], 
  28.                             -1 #  k<0 for diagonals below the main diagonal.
  29.                           )+\
  30.                   np.diag( 1+self.beta[1:self.M],
  31.                           0 # default k=0 for diagonals
  32.                          )+\
  33.                   np.diag( self.gamma[1:self.M-1],
  34.                           1 # k>0 for diagonals above the main diagonal
  35.                          )
  36.     def traverse_grid(self):
  37.         """ Solve using linear systems of equations """
  38.         P,L,U = linalg.lu( self.M1 )
  39.         for j in reversed( range(self.N) ): # 0~t~N-1 ==> N-1~t~0
  40.             # L @ y = B ==> 
  41.             # L @ y = M_2 @ f_j+1 ==> y
  42.             y = linalg.solve( L,  # M_1 @ f_j = M_2 @ f_j+1
  43.                               np.dot( self.M2, self.grid[1:self.M, j+1] )
  44.                             )
  45.             # U @ x = y ==>x
  46.             x = linalg.solve( U,
  47.                               y
  48.                             )
  49.             self.grid[1:self.M, j] = x
  50.             
  51.         print( 'Option price grid', self.grid)

  1. option = FDCnEu(140, 150, r=0.04, T=6/12, 
  2.                       sigma=0.25, Smax=200, M=200, N=2, is_put=False)
  3. print("Call option price:",option.price())

 

  1. option = FDCnEu(140, 150, r=0.04, T=6/12, 
  2.                       sigma=0.25, Smax=200, M=200, N=2, is_put=True)
  3. print("Put option price:",option.price())

 

      Using the same examples that we used with the explicit and implicit methods, we can price a European put option using the Crank-Nicolson method for different time point intervals:

  1. option = FDCnEu(50, 50, r=0.1, T=5./12.,
  2.     sigma=0.4, Smax=100, M=100, N=1000, is_put=True)
  3. print(option.price())

  1. option = FDCnEu(50, 50, r=0.1, T=5./12., 
  2.     sigma=0.4, Smax=100, M=80, N=100, is_put=True)
  3. print(option.price())

 

     From the observed values, the Crank-Nicolson method

  • not only avoids the instability issue we saw in the explicit scheme,
  • but also converges faster than both the explicit and implicit methods.
  • The implicit method requires more iterations, or bigger values of N, to produce values close to those of the Crank-Nicolson method.

Consistency, Stability and Convergence

     When discussing effectiveness of different finite difference methods, we should consider three fundamental properties, which are consistency, stability, and convergence. Before that, we have to notice that though numerical method is good for solving PDE, it also brings the truncation error. Two important sources of error are the truncation error in discretization of stock price and time space.

  • stability: A finite difference scheme is said to be stable if the difference between numerical solution and the true solution will not grow unboundedly不会无限增长​​​​​​​ as dS and dt approach zero. The stability of numerical scheme is achieved if the spectral radius of A less than 1.
  • consistency: A finite difference scheme is consistent if the difference between PDE and finite differential equation vanishes as dS and dt approach zero. That is, the truncation error tends to 0 as the mesh get infinitely finer.
  • convergence: A finite difference scheme is said to be convergent if solutions derived from a finite difference equation converge point-wise to the true solutions of the original PDE逐点收敛到原始 PDE 的真解 as dS and dt approach zero.

Pricing exotic barrier options定价异国障碍期权

     Finite differences are especially useful in pricing exotic options. The nature of the option will dictate(/ ˈdɪkteɪt 规定) the specifications of the boundary conditions.

     In this section, we will take a look at an example of pricing a down-and-out barrier option with the Crank-Nicolson method of finite differences. Due to its relative complexity, other analytical methods, such as Monte Carlo methods, are usually employed in favor of finite difference schemes.

A down-and-out option一种障碍期权

     Let's take a look at an example of a down-and-out option. At any time during the life of the option, should the underlying asset price fall below an S barrier barrier price, the option is considered worthless价格低于某个障碍价格时就自动到期无效. Since, in the grid, the finite difference scheme represents all the possible price points, we only need to consider nodes with the following price range:

We can then set up the boundary conditions as follows: 

pricing down-and-out-options using the Crank-Nicolson method of finite differences

     Let's create a class named FDCnDo that inherits from the FDCnEu class we discussed earlier. We will take into account the barrier price in the constructor method, while leaving the rest of the Crank-Nicolson implementation in the FDCnEu class unchanged:

  1. import numpy as np
  2.  
  3. # Price a down-and-out option by the 
  4. # Crank-Nicolson method of finite differences.
  5.  
  6. class FDCnDo( FDCnEu ):
  7.     def __init__( self, S0, K, r=0.05, T=1, sigma=0,
  8.                   Sbarrier=0, Smax=1, M=1, N=1, is_put=False
  9.                 ):
  10.         # inherit from the 'FDCnEu' class
  11.         # python 2.x
  12.         super( FDCnDo, self ).__init__( 
  13.                   S0, K, r, T, sigma,
  14.                   Smax, M, N, is_put
  15.                 )
  16.         self.barrier = Sbarrier
  17.         """
  18.         https://numpy.org/doc/stable/reference/generated/numpy.linspace.html
  19.         numpy.linspace( start, stop, num=50, endpoint=True,
  20.                         retstep=False, dtype=None, axis=0 )
  21.         np.linspace(2.0, 3.0, num=5)
  22.         array([2.  , 2.25, 2.5 , 2.75, 3.  ])
  23.         """
  24.         self.boundary_conds = np.linspace( Sbarrier, Smax, M+1 )
  25.         self.i_values = self.boundary_conds/self.dS
  26.         
  27.     @property
  28.     def dS(self):
  29.         return (self.Smax-self.barrier)/float(self.M)

     Let's consider an example of a down-and-out option. The underlying stock price is $50 with a volatility of 40 percent. The strike price of the option is $50 with an expiration time of five months(5/12). The risk-free rate is 10 percent. The barrier price is $40.

     We can price a call option and a put down-and-out option with Smax as 100, M as 120, and N as 500

  1. option = FDCnDo(50, 50, r=0.1, T=5./12., sigma=0.4, 
  2.                 Sbarrier=40, Smax=100, M=120, N=500)
  3. print( 'Call option price:', option.price() )

  1. option = FDCnDo(50, 50, r=0.1, T=5./12., sigma=0.4, 
  2.                 Sbarrier=40, Smax=100, M=120, N=500, is_put=True)
  3. print( 'Put option price:', option.price() )

 

The prices of the down-and-out call and put options are $5.4916 and $0.5414, respectively.

Pricing American options with finite differences

     So far, we have priced European options and exotic options. Due to the probability of an early exercise nature in American options, pricing such options is less straightforward. An iterative procedure is required in the implicit Crank-Nicolson method, where the payoffs from earlier exercises in the current period take into account the payoffs of an earlier exercise in the prior period. The Gauss-Siedel iterative method(https://blog.csdn.net/Linli522362242/article/details/125546725) is proposed in the pricing of American options in the Crank-Nicolson method. 

     Recall that in https://blog.csdn.net/Linli522362242/article/details/125546725, The Importance of Linearity in Finance, we covered the Gauss-Siedel method of solving systems of linear equations in the form of Ax=B. Here, the matrix A is decomposed into A=L+U, where L is a lower triangular matrix and U is an upper triangular matrix. Let's take a look at an example of a 4 x 4 matrix, A:

The solution is then obtained iteratively, as follows: 

We can adapt the Gauss-Siedel method to our Crank-Nicolson implementation as follows:
This equation satisfies the early exercise privilege equation:根据无套利原理,在(S,t)空间中任意一点的期权价格不能低于其内在价值(即立刻行使权益的所得收益).对于一般的美式期权，意味着==>在计算得到\LARGE f_{i,j-1}后，我们可以检测提前行使权益的概率

pricing American options using the Crank-Nicolson method of finite differences 

     Let's create a class named FDCnAm that inherits from the FDCnEu class, which is the Crank-Nicolson method's counterpart for pricing European options. The setup_coefficients method may be reused, while overriding all other methods for the inclusion of payoffs from an earlier exercise, if any. 

     The constructor __init__() and the setup_boundary_conditions() methods are given in the FDCnAm class: 

  1. import numpy as np
  2. import sys
  3.  
  4. """ 
  5. Price an American option by the Crank-Nicolson method 
  6. """
  7.  
  8. class FDCnAM( FDCnEu ):
  9.     def __init__(self, S0, K, r=0.05, T=1, sigma=0, 
  10.                  Smax=1, M=1, N=1, 
  11.                  omega=1, tol=0, is_put=False):
  12.         # python 2.x
  13. #         super( FDCnAm, self ).__init__(S0, K, r=r, T=T, sigma=sigma, 
  14. #                                        Smax=Smax, M=M, N=N, 
  15. #                                        is_put=is_put)
  16.         # python 3.x
  17.         super().__init__(S0, K, r=r, T=T, sigma=sigma, 
  18.                          Smax=Smax, M=M, N=N, 
  19.                          is_put=is_put)
  20.         self.omega = omega
  21.         self.tol = tol
  22.         self.i_values = np.arange( self.M+1 ) # Smin~Smax: price level indices
  23.         self.j_values = np.arange( self.N+1 ) # 0~T: time point indices
  24.         
  25.     def setup_boundary_conditions( self ):
  26.         if self.is_call:
  27.             # self.boundary_conds = np.linspace( 0, Smax, self.M+1 )
  28.             # so self.boundary_conds[0] == 0s
  29.             self.payoffs = np.maximum( 0, self.boundary_conds[1:self.M]-self.K )
  30.         else:
  31.             self.payoffs = np.maximum( 0, self.K - self.boundary_conds[1:self.M] )
  32.             
  33.         self.past_values = self.payoffs
  34.         self.boundary_values = self.K * np.exp( -self.r*self.dt*(self.N-self.j_values) )

 

 

      At the start boundaries where the index is 0,the payoffs are calculated with the alpha values omitted. Implement the calculate_payoff_start_boundary() method inside the class:

  1.     def calculate_payoff_start_boundary( self, rhs, old_values ):
  2. #         self.i_values = np.arange( self.M ) # Smin~Smax: price level indices
  3. #         self.j_values = np.arange( self.N ) # 0~T: time point indices        
  4. #     def setup_coefficients(self):
  5. #         self.alpha = 0.25 * self.dt * ( (self.sigma**2)*(self.i_values**2) -\
  6. #                                         self.r * self.i_values
  7. #                                       )
  8. #         self.beta = -0.5 * self.dt * ( (self.sigma**2)*(self.i_values**2) +\
  9. #                                         self.r
  10. #                                     )
  11. #         self.gamma = 0.25 * self.dt * ( (self.sigma**2)*(self.i_values**2) +\
  12. #                                         self.r * self.i_values
  13. #                                       )
  14.         # S[0]~S[M] stock price levels, Smin~Smax
  15.         # old_value.shape : M-1, old_value[0]~old_value[M-2] : f[1]~f[M-1] : S[1]~S[M-1]
  16.         #                                    rhs[0]~rhs[M-2] : f[1]~f[M-1]
  17.         #                        self.beta[0]~[M-1]
  18.         payoff = old_values[0] + self.omega/(1-self.beta[1]) *\
  19.                                  ( rhs[0] -\
  20.                                    ( 1-self.beta[1] ) * old_values[0] +\
  21.                                    self.gamma[1] * old_values[1]
  22.                                  )
  23.         # self.payoffs : option intrinsic value
  24.         return max( self.payoffs[0], payoff )

     At the end boundary where the last index is, the payoffs are calculated with the gamma values omitted. Implement the calculate_payoff_end_boundary() method inside the class:

  1.     def calculate_payoff_end_boundary( self, rhs, old_values, new_values ):
  2.         # old_value.shape : M-1, 0~M stock price levels
  3.         #           f_M-1           rhs_M-1
  4.         payoff = old_values[-1] + self.omega/(1-self.beta[-2]) *\
  5.                                   ( rhs[-1] +\
  6.                                     self.alpha[-2] * new_values[-2] -\
  7.                                     ( 1-self.beta[-2] ) * old_values[-1]
  8.                                   )
  9.         # self.payoffs : option intrinsic value
  10.         return max( self.payoffs[-1], payoff )

     For payoffs that are not at the boundaries, the payoffs are calculated by taking into account the alpha and gamma values. Implement the calculate_payoff() method inside the class:

  1.     def calculate_payoff( self, i, rhs, old_values, new_values ):
  2.         # i start from 2
  3.         payoff = old_values[i] + self.omega/(1-self.beta[i+1]) *\
  4.                                  ( rhs[i] +\
  5.                                    self.alpha[i+1] * new_values[i-1] -\
  6.                                    ( 1-self.beta[i+1] ) * old_values[i] +\
  7.                                    self.gamma[i+1] * old_values[i+1]
  8.                                  )
  9.         return max( self.payoffs[i], payoff )

      Next, implement the traverse_grid() method in the same class:

     The tolerance parameter is used in the Gauss-Siedel method as the convergence criterion. The omega variable is the over-relaxation parameter过松弛参数. Higher omega values provide faster convergence, but this also comes with higher possibilities of the algorithm not converging.

 https://blog.csdn.net/Linli522362242/article/details/125546725
x or \LARGE f<==<==
self.boundary_values = self.K * np.exp( -self.r*self.dt*(self.N-self.j_values) )

  1.     def traverse_grid( self ):
  2.         """ Solve using linear systems of equations """
  3.         aux = np.zeros( self.M-1 ) # 1~M-1 ==> [0,M-2]
  4.         new_values = np.zeros( self.M-1 )
  5.         
  6.         # Gauss-Siede
  7.         for j in reversed( range(self.N) ): # T~0
  8.             aux[0] = self.alpha[1] * ( self.boundary_values[j] +\
  9.                                        self.boundary_values[j+1] 
  10.                                      )
  11.             rhs = np.dot( self.M2, self.past_values ) + aux
  12.             
  13.             # Gauss-Siede
  14.             old_values = np.copy( self.past_values ) # Gauss-Siede
  15.             error = sys.float_info.max
  16.             
  17.             while error > self.tol:
  18.                 new_values[0] = self.calculate_payoff_start_boundary( rhs, old_values )
  19.                 
  20.                 for i in range( self.M-2 )[1:]: # S_1~S[M-3]
  21.                     new_values[i] = self.calculate_payoff( i, rhs, old_values, new_values )
  22.                 
  23.                 new_values[-1] = self.calculate_payoff_end_boundary( rhs, old_values, new_values )
  24.                 
  25.                 error = np.linalg.norm( new_values - old_values )
  26.                 
  27.                 old_values = np.copy( new_values )
  28.                 self.past_values = np.copy( new_values )
  29.                 
  30.         self.values = np.concatenate(( [ self.boundary_values[0] ],
  31.                                        new_values, 
  32.                                        [0]
  33.                                      ), axis=0)

     In each iterative procedure of the while loop, the payoffs are calculated while taking into account the start and end boundaries. Furthermore, new_values are constantly replaced with new payoff calculations based on existing( new_values[i-1] ) and previous values( old_values ).

     Since the new variable, values , contains our terminal payoff values as a one-dimensional array, override the parent interpolate method to account for this change with the following code: 

  1.     def interpolate( self ):
  2.         # Use linear interpolation on final values as 1D array
  3.         """
  4.             Use piecewise linear interpolation on the initial
  5.             grid column to get the closest price at S0.
  6.             
  7.         numpy.interp(x, xp, fp, left=None, right=None, period=None)
  8.         left: Value to return for x < xp[0], default is fp[0]
  9.         right: Value to return for x > xp[-1], default is fp[-1]
  10.         One-dimensional linear interpolation for monotonically 
  11.         increasing sample points.
  12.         Returns the one-dimensional piecewise linear interpolant to
  13.         a function with given discrete data points (xp, fp), 
  14.         evaluated at x.
  15.         xp = [1, 2, 3]
  16.         fp = [3, 2, 0]
  17.         np.interp(2.5, xp, fp) ==> 1.0
  18.         
  19.         self.boundary_conds = np.linspace( 0, Smax, self.M+1 )
  20.         # grid : M+1 stock price levels X N+1 time points
  21.         """
  22.         # self.boundary_conds = np.linspace( 0, Smax, self.M+1 )
  23.         return np.interp( self.S0,
  24.                           self.boundary_conds, self.values
  25.                         ) # lookup grid to return option present value

Put all together: 

  1. import numpy as np
  2. import sys
  3.  
  4. """ 
  5. Price an American option by the Crank-Nicolson method 
  6. """
  7.  
  8. class FDCnAM( FDCnEu ):
  9.     def __init__(self, S0, K, r=0.05, T=1, sigma=0, 
  10.                  Smax=1, M=1, N=1, 
  11.                  omega=1, tol=0, is_put=False):
  12.         # python 2.x
  13. #         super( FDCnAm, self ).__init__(S0, K, r=r, T=T, sigma=sigma, 
  14. #                                        Smax=Smax, M=M, N=N, 
  15. #                                        is_put=is_put)
  16.         # python 3.x
  17.         super().__init__(S0, K, r=r, T=T, sigma=sigma, 
  18.                          Smax=Smax, M=M, N=N, 
  19.                          is_put=is_put)
  20.         self.omega = omega
  21.         self.tol = tol
  22.         self.i_values = np.arange( self.M+1 ) # Smin~Smax: price level indices
  23.         self.j_values = np.arange( self.N+1 ) # 0~T: time point indices
  24.         
  25.     def setup_boundary_conditions( self ):
  26.         if self.is_call:
  27.             # self.boundary_conds = np.linspace( 0, Smax, self.M+1 )
  28.             # so self.boundary_conds[0] == 0s
  29.             self.payoffs = np.maximum( 0, self.boundary_conds[1:self.M]-self.K )
  30.         else:
  31.             self.payoffs = np.maximum( 0, self.K - self.boundary_conds[1:self.M] )
  32.             
  33.         self.past_values = self.payoffs
  34.         self.boundary_values = self.K * np.exp( -self.r*self.dt*(self.N-self.j_values) )
  35.         
  36.     def calculate_payoff_start_boundary( self, rhs, old_values ):
  37. #         self.i_values = np.arange( self.M ) # Smin~Smax: price level indices
  38. #         self.j_values = np.arange( self.N ) # 0~T: time point indices        
  39. #     def setup_coefficients(self):
  40. #         self.alpha = 0.25 * self.dt * ( (self.sigma**2)*(self.i_values**2) -\
  41. #                                         self.r * self.i_values
  42. #                                       )
  43. #         self.beta = -0.5 * self.dt * ( (self.sigma**2)*(self.i_values**2) +\
  44. #                                         self.r
  45. #                                     )
  46. #         self.gamma = 0.25 * self.dt * ( (self.sigma**2)*(self.i_values**2) +\
  47. #                                         self.r * self.i_values
  48. #                                       )
  49.         # S[0]~S[M] stock price levels, Smin~Smax
  50.         # old_value.shape : M-1, old_value[0]~old_value[M-2] : f[1]~f[M-1] : S[1]~S[M-1]
  51.         #                                    rhs[0]~rhs[M-2] : f[1]~f[M-1]
  52.         #                        self.beta[0]~[M-1]
  53.         payoff = old_values[0] + self.omega/(1-self.beta[1]) *\
  54.                                  ( rhs[0] -\
  55.                                    ( 1-self.beta[1] ) * old_values[0] +\
  56.                                    self.gamma[1] * old_values[1]
  57.                                  )
  58.         # self.payoffs : option intrinsic value
  59.         return max( self.payoffs[0], payoff )
  60.     
  61.     def calculate_payoff_end_boundary( self, rhs, old_values, new_values ):
  62.         # old_value.shape : M-1, 0~M stock price levels
  63.         #           f_M-1           rhs_M-1
  64.         payoff = old_values[-1] + self.omega/(1-self.beta[-2]) *\
  65.                                   ( rhs[-1] +\
  66.                                     self.alpha[-2] * new_values[-2] -\
  67.                                     ( 1-self.beta[-2] ) * old_values[-1]
  68.                                   )
  69.         # self.payoffs : option intrinsic value
  70.         return max( self.payoffs[-1], payoff )
  71.     
  72.     def calculate_payoff( self, i, rhs, old_values, new_values ):
  73.         # i start from 2
  74.         payoff = old_values[i] + self.omega/(1-self.beta[i+1]) *\
  75.                                  ( rhs[i] +\
  76.                                    self.alpha[i+1] * new_values[i-1] -\
  77.                                    ( 1-self.beta[i+1] ) * old_values[i] +\
  78.                                    self.gamma[i+1] * old_values[i+1]
  79.                                  )
  80.         return max( self.payoffs[i], payoff )
  81.     
  82.     def traverse_grid( self ):
  83.         """ Solve using linear systems of equations """
  84.         aux = np.zeros( self.M-1 ) # 1~M-1 ==> [0,M-2]
  85.         new_values = np.zeros( self.M-1 )
  86.         
  87.         # Gauss-Siede
  88.         for j in reversed( range(self.N) ): # T~0
  89.             aux[0] = self.alpha[1] * ( self.boundary_values[j] +\
  90.                                        self.boundary_values[j+1] 
  91.                                      )
  92.             rhs = np.dot( self.M2, self.past_values ) + aux
  93.             
  94.             # Gauss-Siede
  95.             old_values = np.copy( self.past_values ) # Gauss-Siede
  96.             error = sys.float_info.max
  97.             
  98.             while error > self.tol:
  99.                 new_values[0] = self.calculate_payoff_start_boundary( rhs, old_values )
  100.                 
  101.                 for i in range( self.M-2 )[1:]: # S_1~S[M-3]
  102.                     new_values[i] = self.calculate_payoff( i, rhs, old_values, new_values )
  103.                 
  104.                 new_values[-1] = self.calculate_payoff_end_boundary( rhs, old_values, new_values )
  105.                 
  106.                 error = np.linalg.norm( new_values - old_values )
  107.                 
  108.                 old_values = np.copy( new_values )
  109.                 self.past_values = np.copy( new_values )
  110.                 
  111.         self.values = np.concatenate(( [ self.boundary_values[0] ],
  112.                                        new_values, 
  113.                                        [0]
  114.                                      ), axis=0)
  115.         
  116.     def interpolate( self ):
  117.         # Use linear interpolation on final values as 1D array
  118.         """
  119.             Use piecewise linear interpolation on the initial
  120.             grid column to get the closest price at S0.
  121.             
  122.         numpy.interp(x, xp, fp, left=None, right=None, period=None)
  123.         left: Value to return for x < xp[0], default is fp[0]
  124.         right: Value to return for x > xp[-1], default is fp[-1]
  125.         One-dimensional linear interpolation for monotonically 
  126.         increasing sample points.
  127.         Returns the one-dimensional piecewise linear interpolant to
  128.         a function with given discrete data points (xp, fp), 
  129.         evaluated at x.
  130.         xp = [1, 2, 3]
  131.         fp = [3, 2, 0]
  132.         np.interp(2.5, xp, fp) ==> 1.0
  133.         
  134.         self.boundary_conds = np.linspace( 0, Smax, self.M+1 )
  135.         # grid : M+1 stock price levels X N+1 time points
  136.         """
  137.         # self.boundary_conds = np.linspace( 0, Smax, self.M+1 )
  138.         return np.interp( self.S0,
  139.                           self.boundary_conds, self.values
  140.                         ) # lookup grid to return option present value

     Let's price an American call-and-put option with an underlying asset price of 50 and volatility of 40 percent, a strike price of 50, a risk-free rate of 10 percent, and an expiration date of five months(5./12). We choose a Smax value of 100 , M as 100 , N as 42 , an omega parameter
value of 1.2 , and a tolerance value of 0.001:

  1. am_call_option = FDCnAM( 50, 50, r=0.1, T=5./12., 
  2.                          sigma=0.4, Smax=100, M=100, N=42,
  3.                          omega=1.2, tol=0.001, is_put=False
  4.                        )
  5. print( "American call option price:",am_call_option.price() )

  1. am_put_option = FDCnAM( 50, 50, r=0.1, T=5./12., 
  2.                          sigma=0.4, Smax=100, M=100, N=42,
  3.                          omega=1.2, tol=0.001, is_put=True
  4.                        )
  5. print( "American put option price:",am_put_option.price() )

 

      The prices of the call and put American stock options by using the Crank-Nicolson method are $6.109 and $4.2778, respectively.

  1. am_call_option = FDCnAM( 140, 150, r=0.04, T=6./12., 
  2.                          sigma=0.25, Smax=200, M=200, N=2,
  3.                          omega=1.2, tol=0.001, is_put=False
  4.                        )
  5. print( "American call option price:",am_call_option.price() )

  1. am_put_option = FDCnAM( 140, 150, r=0.04, T=6./12., 
  2.                          sigma=0.25, Smax=200, M=200, N=2,
  3.                          omega=1.2, tol=0.001, is_put=True
  4.                        )
  5. print( "American put option price:",am_put_option.price() )

 
     Since American options can be exercised at any time and European options can only be exercised at maturity, this added flexibility of American options increases their value over European options in certain circumstances.

Putting it all together – implied volatility modeling

     In the option pricing methods we have learned so far, a number of parameters are assumed to be constant许多参数被假定为常数: interest rates, strike prices, dividends, and volatility. Here, the parameter of interest is volatility. In quantitative research, the volatility ratio is used to forecast price trends.

     To derive implied volatilities, we need to refer to https://blog.csdn.net/Linli522362242/article/details/125662545, Nonlinearity in Finance, where we discussed the root-finding methods of nonlinear functions. We will use the bisection method二分法 of numerical procedures in our next example to create an implied volatility curve.

Implied volatilities of the AAPL American put option

     Let's consider the option data of the stock Apple (AAPL), which was gathered at the end of the day on October 3, 2014. These details are provided in the following table. The option expires on December 20, 2014. The prices listed are the mid-points of the bid and ask prices:

     The last traded price of AAPL was 99.62, with an interest rate of 2.48 percent and a dividend yield of 1.82 percent. The American options expire in 78 days.

     Using this information, let's create a new class named ImpliedVolatilityModel that accepts the stock option's parameters in the constructor. If required, import the BinomialLROption class that we created for the LR binomial tree we covered in the earlier section of this chapter, A class for the LR binomial tree option pricing model.

  1. import math
  2. # Stores common attributes of stock option
  3. class StockOption(object):
  4.     def __init__(
  5.         self, S0, K, r =0.05, T=1, N=2, pu=0, pd=0,
  6.         div=0, sigma=0, is_put = False, is_am=False
  7.     ):
  8.         """
  9.         Initialize the stock option base class.
  10.         Defaults to European call unless specified.
  11.         : param S0: initial stock price
  12.         : param K: strike price
  13.         : param r: risk-free interest rate
  14.         : param T: time to maturity
  15.         : param N: number of time steps
  16.         : param pu: probability at up state
  17.         : param pd: probability at down state
  18.         : param div: Dividend yield
  19.         : param is_put: True for a put option,
  20.                         False for a call option
  21.         : param is_am: True for an American option,
  22.                        False for a European option
  23.         """
  24.         self.S0=S0
  25.         self.K=K
  26.         self.r=r
  27.         self.T=T
  28.         self.N = max(1,N)
  29.         self.STs = [] # Declare the stock prices tree
  30.         
  31.         # Optional parameters used by derived classes
  32.         self.pu, self.pd = pu, pd
  33.         self.div = div
  34.         self.sigma = sigma
  35.         self.is_call = not is_put
  36.         self.is_european = not is_am
  37.         # @property:https://www.programiz.com/python-programming/property
  38.     @property
  39.     def dt(self):
  40.         # Single time step, in years
  41.         return self.T/float(self.N)
  42.     
  43.     @property
  44.     def df(self):
  45.         # The discount factor
  46.         return math.exp( -(self.r - self.div)*self.dt )
  47.  
  48. ############################################
  49. import math
  50. import numpy as np
  51.  
  52. # Price a European or American option by the binomial tree
  53. class BinomialTreeOption(StockOption):
  54.     def setup_parameters( self ):
  55.         self.u = 1 + self.pu # Expected value % in the up state, e.g 1+0.2=1.2
  56.         self.d = 1 - self.pd # Expected value % in the down state e.g 1-0.2=0.8
  57.         
  58.         # qu : The risk-neutral (up state) probability 
  59.         self.qu = ( math.exp( (self.r-self.div)*self.dt )-self.d 
  60.                   )/(self.u-self.d)
  61.         # qd : The risk-neutral (down state) probability
  62.         self.qd = 1-self.qu
  63.         
  64.     def init_stock_price_tree(self):
  65.         # Initialize a 2D tree at T=0
  66.         # to store the expected returns of the stock prices for all time steps
  67.         self.STs = [ np.array([self.S0]) ] # [array([50])] # root
  68.         
  69.         # Simulate the possible stock prices path
  70.         # calculate the payoff values 
  71.         # from exercising the option at each period
  72.         for i in range( self.N ): # N: number of time steps(branch_nodes)
  73.             prev_branches = self.STs[-1]     # at last time step
  74.             st = np.concatenate(             # note, cross node: Sud = Sdu
  75.                 ( prev_branches*self.u,      # go up
  76.                   [prev_branches[-1]*self.d] # the last one goes down
  77.                 )# array([60., 40.])
  78.             )    # array([72., 48., 32.])
  79.             
  80.             self.STs.append(st) # Add nodes at each time step
  81.     
  82.     # the intrinsic values of the option at maturity
  83.     def init_payoffs_tree(self):
  84.         if self.is_call:# call option
  85.             return np.maximum( 0, self.STs[self.N]-self.K )
  86.         else:# put option
  87.             return np.maximum( 0, self.K-self.STs[self.N] )
  88.              # e.g. payoffs: [ 0.  4. 20.]
  89.     
  90.     # Returns the maximum payoff values between 
  91.     # exercising the American option early 
  92.     # and not exercising the option at all
  93.     def check_early_exercise( self, payoffs, node ):
  94.         if self.is_call:
  95.             return np.maximum( payoffs, self.STs[node]-self.K )
  96.         else:
  97.             # payoffs: [1.41475309 9.46393007]
  98.             # payoffs at maturity: [-8. 12.]
  99.             # np.maximum ==>[1.41475309 12.]
  100.             return np.maximum( payoffs, self.K-self.STs[node] )
  101.     
  102.     def traverse_tree( self, payoffs ):
  103.         """
  104.         Starting from the time the option expires, traverse
  105.         backwards and calculate discounted payoffs at each node
  106.         """
  107.         for i in reversed( range(self.N) ):
  108.             # The payoffs from NOT exercising the option
  109.             payoffs = ( payoffs[:-1]*self.qu +
  110.                         payoffs[1:]*self.qd
  111.                       ) * self.df # df: The discount factor
  112.                       # * math.exp( -(self.r - self.div)*self.dt )
  113.             print("At the time step",i,", payoffs:",payoffs)
  114.             # [1.41475309 9.46393007]
  115.             # [4.19265428]
  116.             # Payoffs from exercising, for American options
  117.             if not self.is_european:
  118.                 payoffs = self.check_early_exercise( payoffs, i )
  119.                 print("At the time step",i,", max payoffs:",payoffs)
  120.                 
  121.         return payoffs
  122.     
  123.     def begin_tree_traversal( self ):
  124.         # the intrinsic values of the option at maturity
  125.         payoffs = self.init_payoffs_tree() # e.g. payoffs: [ 0. 4. 20.]
  126.         print("the payoffs at maturity", payoffs)
  127.         return self.traverse_tree( payoffs )
  128.     
  129.     ###### Entry point of the pricing implementation ######
  130.     def price(self):
  131.         self.setup_parameters()
  132.         self.init_stock_price_tree()
  133.         # self.STs : [array([50]), array([60., 40.]), array([72., 48., 32.])]
  134.         print("Nodes:",self.STs)
  135.         payoffs = self.begin_tree_traversal()
  136.         #print(payoffs)
  137.         return payoffs[0]
  138.  
  139. ############################################
  140. import math
  141.  
  142. # Price an option by the Leisen-Reimer tree
  143.  
  144. class BinomialLROption( BinomialTreeOption ):
  145.     def pp_2_inversion( self, z, n ):
  146.         # copysign : return the value of the first parameter and the sign of the second parameter
  147.         # math.copysign(4, -1)     ==> -4.0
  148.         # math.copysign(-8, 97.21) ==> 8.0
  149.         # math.copysign(-43, -76)  ==> -43.0
  150.         return 0.5 + math.copysign(1,z) *\
  151.             math.sqrt( 0.25-0.25*math.exp( 
  152.                         -( ( z/( n+1./3.+0.1/(n+1) ) )**2.
  153.                          ) * ( n+1./6. ) )
  154.                      )
  155.     
  156.     def setup_parameters( self ):
  157.         odd_N = self.N if (self.N%2==0) else (self.N+1)
  158.         d1 = ( math.log(self.S0/self.K) + 
  159.               ( (self.r-self.div) + (self.sigma**2)/2. )*self.T
  160.              ) / ( self.sigma * math.sqrt(self.T) )
  161.  
  162. #         d2 = ( math.log(self.S0/self.K) + 
  163. #               ( (self.r-self.div) - (self.sigma**2)/2. )*self.T
  164. #              ) / ( self.sigma * math.sqrt(self.T) )
  165.         d2 = d1-self.sigma*math.sqrt(self.T)
  166.         
  167.         pbar = self.pp_2_inversion( d1, odd_N )
  168.         self.p = self.pp_2_inversion( d2, odd_N )
  169.         
  170.         # self.dt : self.T/float(self.N)
  171.         # self.df : math.exp( -(self.r - self.div)*self.dt )
  172.         # self.u : Expected value in the up state, e.g 1+0.2=1.2
  173.         self.u = 1/self.df * pbar/self.p
  174.         # self.d : Expected value in the down state e.g 1-0.2=0.8
  175.         self.d = (1/self.df - self.p*self.u)/(1-self.p)
  176.         
  177.         # qu : The risk-neutral (up state) probability
  178.         self.qu = self.p
  179.         # qd : The risk-neutral (down state) probability
  180.         self.qd = 1-self.p

     The bisection function we covered in https://blog.csdn.net/Linli522362242/article/details/125662545, Nonlinearity in Finance, is also required.

Import the bisection function we discussed in the previous chapter:

  1. # The bisection method
  2. import numpy as np
  3.  
  4. def bisection(func, a, b, tol=0.1, maxiter=20):
  5.     """
  6.     :param func: The function to solve
  7.     :param a: The x-axis value where f(a)<0
  8.     :param b: The x-axis value where f(b)>0
  9.     :param tol: The precision of the solution
  10.     :param maxiter: Maximum number of iterations
  11.     :return: 
  12.         The x-axis value of the root,
  13.         number of iterations used
  14.     """
  15.     c = (a+b)*0.5 # Declare c as the midpoint [a,b]
  16.     n = 1 # start with 1 iteration
  17.     while n<= maxiter:
  18.         c = (a+b)*0.5
  19.         # Root is found or is very close
  20.         if func(c) == 0 or abs(a-b)*0.5 < tol:
  21.             return c,n
  22.         n += 1
  23.         
  24.         # Decide the side to repeat the steps
  25.         if func(c) < 0:
  26.             a=c
  27.         else:
  28.             b=c
  29.             
  30.     return c,n

     The option_valuation() method accepts the K strike price and the sigma volatility value to compute the value of the option. In this example, we are using the BinomialLROption pricing method.

     The get_implied_volatilities() method accepts a list of strike and option prices to compute the implied volatilities by the bisection method for every price available. Therefore, the length of the two lists must be the same.

The Python code for the ImpliedVolatilityModel class is given as follows:

  1. """
  2. Get implied volatilities from a Leisen-Reimer binomial tree 
  3. using the bisection method as the numerical procedure.
  4. """
  5.  
  6. class ImpliedVolatilityModel( object ):
  7.     def __init__( self, S0, r=0.05, T=1, div=0,
  8.                   N=1, is_put=False ):
  9.         self.S0 = S0
  10.         self.r = r
  11.         self.T = T
  12.         self.div = div
  13.         self.N = N
  14.         self.is_put = is_put
  15.         
  16.     def option_valuation(self, K, sigma):
  17.         # accepts the K strike price and the sigma volatility value 
  18.         # to compute the value of the option.
  19.         """ Use the binomial Leisen-Reimer tree """
  20.         lr_option = BinomialLROption( self.S0, K, r=self.r, T=self.T, N=self.N,
  21.                                       sigma=sigma,
  22.                                       div=self.div, is_put=self.is_put 
  23.                                     )
  24.         return lr_option.price()
  25.     
  26.     def get_implied_volatilities( self, Ks, opt_prices ):
  27.         impvols = []
  28.         
  29.         for i in range( len(strikes) ):
  30.             # Bind f(sigma) for use by the bisection method
  31.             f = lambda sigma: self.option_valuation( Ks[i], sigma ) - opt_prices[i]
  32.             
  33.             impv = bisection(f, 0.01, 0.99, 0.0001, 100)[0]
  34.             # solution : self.option_valuation( Ks[i], sigma ) == opt_prices[i] ==> sigma
  35.             impvols.append(impv)
  36.             
  37.         return impvols

     Using this model, let's find out the implied volatilities of the American put options using this particular set of data:

  1. strikes = [75, 80, 85, 90, 92.5, 95, 97.5, 
  2.            100, 105, 110, 115, 120, 125]
  3. put_prices = [0.16, 0.32, 0.6, 1.22, 1.77, 2.54, 3.55, 
  4.               4.8, 7.75, 11.8, 15.96, 20.75, 25.81]
  5.  
  6. model = ImpliedVolatilityModel( 99.62, r=0.0248, T=78/365., div=0.0182, 
  7.                                 N=77, is_put=True
  8.                               )
  9. impvols_put = model.get_implied_volatilities(strikes, put_prices)
  10. impvols_put


     The implied volatility values are now stored in the impvols_put variable as a list object. Let's plot these values against the strike prices to obtain an implied volatility curve:

  1. %matplotlib inline
  2. import matplotlib.pyplot as plt
  3.  
  4. plt.plot( strikes, impvols_put )
  5. plt.xlabel('Strike Prices')
  6. plt.ylabel('Implied Volatilities')
  7. plt.title('AAPL Put Implied Volatilities expiring in 78 days')
  8. plt.show()

     This would give us the volatility smile, as shown in the following diagram. Here, we have modeled an LR tree with 77 steps, with each step representing one day: 

 
     Of course, pricing an option daily may not be ideal since markets change by fractions of a millisecond. We used the bisection method to solve the implied volatility as implied by the binomial tree, as opposed to the realized volatility values directly observed from market prices.

     Should we fit this curve against a polynomial curve to identify potential arbitrage opportunities? Or extrapolate the curve to derive further insights on potential opportunities from implied volatilities of far out-of-the-money
(There are also options traded that are far out-of-the-money (index level (PRICE) much lower than option strike(STRIKE)).​​​​​​​)
and in-the-money
(there are call options traded and quoted that are far in-the-money (index level (PRICE) much higher than option strike (STRIKE)).)
options? Well, these questions are for option traders like yourself to find out! 

Summary

     In this chapter, we looked at a number of numerical procedures in derivative pricing, the most common being options. One such procedure is the use of trees, with binomial trees being the simplest structure to model asset information, where one node extends to two other nodes in each time step, representing an up state and a down state, respectively. In trinomial trees, each node extends to three other nodes in each time step, representing an up state, a down state, and a state with no movement, respectively. As the tree traverses upwards, the underlying asset is computed and represented at each node. The option then takes on the structure of this tree and, starting from the terminal payoffs, the tree traverses backward and toward the root, which converges to the current discounted option price. Besides binomial and trinomial trees, trees can take on the form of the CRR, Jarrow-Rudd, Tian, or LR parameters.

     By adding another layer of nodes around our tree, we introduced additional information from which we can derive the Greeks, such as the delta and gamma, without incurring additional computational costs.

     Lattices were introduced as a way of saving storage costs over binomial and trinomial trees. In lattice pricing, nodes with new information are saved only once and reused later on nodes that require no change in the information.

     We also discussed the finite difference schemes in option pricing, consisting of terminal and boundary conditions. From the terminal conditions, the grid traverses backward in time using the explicit method, implicit method, and the Crank-Nicolson method. Besides pricing European and American options, finite difference pricing schemes can be used to price exotic options, where we looked at an example of pricing a down-and-out barrier option.

     By importing the bisection root-finding method learned about in https://blog.csdn.net/Linli522362242/article/details/125662545, Nonlinearity in Finance, and the binomial LR tree model in this chapter, we used market prices of an American option to create an implied volatility curve for further studies.

In the next chapter, we will take a look at modeling interest rates and derivatives.

声明：本文内容由网友自发贡献，转载请注明出处：【wpsshop博客】
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号