Dividend models for BS formula.

It would be quite handy to know what is the true price of a european option given cash dividends. Without dividends we have the Black and Scholes analytical formula. With dividends, we can have a lower and a upper bound.

For the lower bound, we use the approach I called "Hull":
The stock $S$ is decomposed into a risky component $S^*$ and a certain component $D$ which is the present value of all dividends to be paid during the life of the option.
Since they will be paid before the option expired, the method prices the option using only $S^*$.
We argued that it was underestimating the price because the volatility is applied to $S^*$ until the dividends are paid and it can be quite smaller than $S$. The sooner the dividends are paid, the smaller the error though. Indeed if the dividends would be paid right now, it would be correct to price the option with $S^*$

For the upper bound, one can use the forward method:
If we had the dividend paid right at expiration of the option, payoff of the call would read 'max ((S-D)-K,0)' and the one of the put 'max((-(S-D)+K,0)'. In that case we could price the option using S and using a strike K+D instead of K.
The forward method extended this: all dividends are grown until expiration (if time of the dividend is $t_1$ and expiration $T$ we use $e^{r(T-t_1)}$ and added to the strike.
The volatility is here applied to the full stock $S$ even after dividend are paid. So it overestimate the price of an option.

There is some literature on how to correct the volatility we use with the Hull or the forward methods to get a good price but as soon as you have several dividends it becomes messy.

Find below the code for both adaptation of the BS formula for dividend options using the Hull and the forward method:

1. def BlackScholes2b(option):  
2.     [S0,K,sigma,T,r,divi,american]=[float(option.S0),float(option.K),option.sigma,float(option.T),option.r,option.divi,option.american]  
3.     #Hull method, dividend are discounted and removed from the spot  
4.     
5.     def phi(z):  
6.         return np.exp(-.5 * z * z) / (np.sqrt( 2 * np.pi ))  
7.       
8.     dividends=[[],[]]  
9.     pvdividends=0  
10.       
11.     if (np.size(divi)>0 and divi[0][0]<T) :  
12.         lastdiv=np.nonzero(np.array(divi[0][:])<=T)[0][-1]  
13.         dividends[0]=divi[0][:lastdiv+1]          
14.         dividends[1]=divi[1][:lastdiv+1]   
15.       
16.     if np.size(dividends)>0:  
17.         for index,dv in enumerate(list(dividends[1])):  
18.             pvdividends=pvdividends+dv*np.exp(-r*dividends[0][index])  
19.       
20.     S0=S0-pvdividends  
21.     td1=d1(S0, K, r, sigma, T)  
22.     td2=d2(S0, K, r, sigma, T)  
23.       
24.     option.call.price=S0 * ss.norm.cdf(td1) - K * np.exp(-r * T) * ss.norm.cdf(td2)  
25.     option.put.price=K * np.exp(-r * T) * ss.norm.cdf(-td2) - S0 * ss.norm.cdf(-td1)  
26.   
27.     option.call.delta=ss.norm.cdf(td1)  
28.     option.put.delta=ss.norm.cdf(td1)-1  
29.       
30.     NPrime=phi(td1)  
31.     option.call.gamma=NPrime/(S0*sigma*np.sqrt(T))  
32.     option.put.gamma=option.call.gamma  
33.       
34.     option.call.theta=(NPrime)*(-S0*sigma*0.5/np.sqrt(T))-r*K * np.exp(-r * T) * ss.norm.cdf(td2)  
35.     option.put.theta=(NPrime)*(-S0*sigma*0.5/np.sqrt(T))+r*K * np.exp(-r * T) * ss.norm.cdf(-td2)  
36.           
37.     return option  
38.       
39.       
40. def BlackScholes2c(option):  
41.     [S0,K,sigma,T,r,divi,american]=[float(option.S0),float(option.K),option.sigma,float(option.T),option.r,option.divi,option.american]  
42.     #Forward method, dividends are grown and added to the strike  
43.     
44.     def phi(z):  
45.         return np.exp(-.5 * z * z) / (np.sqrt( 2 * np.pi ))  
46.       
47.     dividends=[[],[]]  
48.     pvdividends=0  
49.       
50.       
51.     if (np.size(divi)>0 and divi[0][0]<T) :  
52.         lastdiv=np.nonzero(np.array(divi[0][:])<=T)[0][-1]  
53.         dividends[0]=divi[0][:lastdiv+1]          
54.         dividends[1]=divi[1][:lastdiv+1]   
55.       
56.     if np.size(dividends)>0:  
57.        for index,dv in enumerate(list(dividends[1])):  
58.             pvdividends=pvdividends+dv*np.exp(r*(T-dividends[0][index]))  
59.     
60.     K=K+pvdividends  
61.           
62.     td1=d1(S0, K, r, sigma, T)  
63.     td2=d2(S0, K, r, sigma, T)  
64.       
65.     option.call.price=S0 * ss.norm.cdf(td1) - K * np.exp(-r * T) * ss.norm.cdf(td2)  
66.     option.put.price=K * np.exp(-r * T) * ss.norm.cdf(-td2) - S0 * ss.norm.cdf(-td1)  
67.   
68.     option.call.delta=ss.norm.cdf(td1)  
69.     option.put.delta=ss.norm.cdf(td1)-1  
70.       
71.     NPrime=phi(td1)  
72.     option.call.gamma=NPrime/(S0*sigma*np.sqrt(T))  
73.     option.put.gamma=option.call.gamma  
74.       
75.     option.call.theta=(NPrime)*(-S0*sigma*0.5/np.sqrt(T))-r*K * np.exp(-r * T) * ss.norm.cdf(td2)  
76.     option.put.theta=(NPrime)*(-S0*sigma*0.5/np.sqrt(T))+r*K * np.exp(-r * T) * ss.norm.cdf(-td2)  
77.           
78.     return option