-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathv1_option_pricing_cy.pyx
178 lines (168 loc) · 7.72 KB
/
v1_option_pricing_cy.pyx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
import numpy as np
cimport numpy as np
cimport cython
cdef extern from "math.h":
double exp(double x)
double pow(double base, double exponent)
double sqrt(double x)
double fmax(double x, double y)
# Disable bounds checking for performance
@cython.boundscheck(False)
# Function to calculate option price using binomial model
def option_binomial(
bint model, # 1 American, 0 European
float flag, # Flag to indicate whether it's a call or put option
float S, # Initial stock price
float X, # Strike price
float r, # Risk-free rate
float sigma, # Volatility
float t, # Time to expiration
int steps, # Number of steps in the binomial tree
float div_yield # Dividend yield
):
"""
This function calculates the option price using a binomial model.
Parameters:
american
flag (float): Flag to indicate whether it's a call or put option
S (float): Initial stock price
X (float): Strike price
r (float): Risk-free rate
sigma (float): Volatility
t (float): Time to expiration
steps (int): Number of steps in the binomial tree
div_yield (float): Dividend yield
Returns:
float: The calculated option price
"""
# Define local variables
cdef int cstep
cdef int x
cdef int step
cdef int i
cdef float R = exp((r - div_yield) * (t/steps)) # Discount factor per step, adjusted for dividend yield
cdef float R_no_div_yield = exp(r*(t/steps))
cdef float Rinv = 1.0/R_no_div_yield # Inverse of discount factor
cdef float u = exp(sigma * sqrt(t / steps)) # Upward movement factor
cdef float uu = u * u # Square of upward movement factor
cdef float d = 1.0/u # Downward movement factor
cdef float p_up = (R - d) / (u - d) # Probability of upward movement
cdef float p_down = 1-p_up # Probability of downward movement
cdef np.ndarray[np.double_t, ndim=1] prices # Array to store stock prices
cdef np.ndarray[np.double_t, ndim=1] option_values # Array to store option values
prices = np.zeros(steps + 1, dtype=np.double) # Initialize prices array
option_values = np.zeros(steps + 1, dtype=np.double) # Initialize option values array
prices[0] = S * pow(d, steps) # Calculate initial stock price
for i in range(1, steps + 1):
prices[i] = uu * prices[i-1] # Calculate stock price for each step
for i in range(steps+1):
option_values[i] = max(0., flag * (prices[i]-X)) # Calculate option value for each step
for step in range(steps-1, -1, -1):
cstep = step
for i in range(cstep+1):
# Update option value based on binomial model
option_values[i] = (p_up * option_values[i+1] + p_down * option_values[i])*Rinv
prices[i] = d * prices[i+1] # Update stock price
# Update option value based on exercise decision
if model:
option_values[i] = max(option_values[i], flag*(prices[i]-X))
return option_values[0] # Return the option price
# Disable bounds checking for performance
@cython.boundscheck(False)
# Function to calculate option price with discrete dividends using binomial model
def discrete_divs_cy(
bint model, # 1 American, 2 European
float flag, # Flag to indicate whether it's a call or put option
float S, # Initial stock price
float X, # Strike price
float r, # Risk-free rate
float sigma, # Volatility
float t, # Time to expiration
int steps, # Number of steps in the binomial tree
np.ndarray[np.double_t, ndim=1] div_times, # Array of dividend times
np.ndarray[np.double_t, ndim=1] div_amts, # Array of dividend amounts
float div_yield # Dividend yield
):
"""
This function calculates the option price with discrete dividends using a binomial model.
Parameters:
model (bint): Model is American or European
flag (float): Flag to indicate whether it's a call or put option
S (float): Initial stock price
X (float): Strike price
r (float): Risk-free rate
sigma (float): Volatility
t (float): Time to expiration
steps (int): Number of steps in the binomial tree
div_times (np.ndarray): Array of dividend times
div_amts (np.ndarray): Array of dividend amounts
div_yield (float): Dividend yield
Returns:
float: The calculated option price
"""
# Define local variables
cdef int n_dividends = div_times.shape[0]
if n_dividends == 0:
# If no dividends, use simple binomial model
return option_binomial(model, flag, S, X, r, sigma, t, steps, div_yield)
cdef int steps_before = <int> (steps*(div_times[0]/t))
#print("steps_before:", steps_before)
if steps_before < 0:
steps_before = 0
if steps_before > steps:
steps_before = steps-1
cdef double value_alive
cdef int cstep
cdef int x
cdef int step
cdef int i
cdef float R = exp((r - div_yield) * (t/steps)) # Discount factor per step, adjusted for dividend yield
cdef float R_no_div_yield = exp(r*(t/steps))
cdef float Rinv = 1.0/R_no_div_yield # Inverse of discount factor
cdef float u = exp(sigma * sqrt(t/steps)) # Upward movement factor
cdef float uu = u * u # Square of upward movement factor
cdef float d = 1.0/u # Downward movement factor
cdef float p_up = (R-d)/(u-d) # Probability of upward movement
cdef float p_down = 1-p_up # Probability of downward movement
cdef double dividend_amount = div_amts[0]
cdef np.ndarray[np.double_t, ndim=1] tmp_dividend_times
cdef np.ndarray[np.double_t, ndim=1] tmp_dividend_amts
cdef np.ndarray[np.double_t, ndim=1] prices # Array to store stock prices
cdef np.ndarray[np.double_t, ndim=1] option_values # Array to store option values
if n_dividends > 1:
tmp_dividend_times = np.zeros(n_dividends-1, dtype=np.double)
tmp_dividend_amts = np.zeros(n_dividends-1, dtype=np.double)
for i in range(0, n_dividends-1, 1):
tmp_dividend_times[i] = div_times[i-1]
tmp_dividend_amts[i] = div_amts[i-1]
prices = np.zeros(steps_before+1, dtype=np.double)
option_values = np.zeros(steps_before+1, dtype=np.double)
prices[0]=S*pow(d, steps_before)
for i in range(1, steps_before+1):
prices[i] = uu * prices[i-1]
for i in range(steps_before+1):
value_alive = discrete_divs_cy(model, flag, prices[i]-dividend_amount, X, r, sigma, t-div_times[0], steps-steps_before, tmp_dividend_times, tmp_dividend_amts, div_yield)
option_values[i] = max(value_alive, flag * (prices[i]-X))
for step in range(steps_before-1, -1, -1):
cstep = step
for i in range(cstep+1):
option_values[i] = (p_up * option_values[i+1] + p_down * option_values[i])*Rinv
prices[i] = d * prices[i+1]
if model:
option_values[i] = max(option_values[i], flag*(prices[i]-X))
else:
prices = np.zeros(steps_before+1, dtype=np.double)
option_values = np.zeros(steps_before+1, dtype=np.double)
prices[0]=S*pow(d, steps_before)
for i in range(1, steps_before+1):
prices[i] = uu * prices[i-1]
for i in range(steps_before+1):
value_alive = option_binomial(model, flag, prices[i]-dividend_amount, X, r, sigma, t-div_times[0], steps-steps_before, div_yield)
option_values[i] = max(value_alive, flag * (prices[i]-X))
for step in range(steps_before-1, -1, -1):
cstep = step
for i in range(cstep+1):
option_values[i] = (p_up * option_values[i+1] + p_down * option_values[i])*Rinv
prices[i] = d * prices[i+1]
option_values[i] = max(option_values[i], flag*(prices[i]-X))
return option_values[0] # Return the option price