Coverage for rivapy/models/heston.py: 96%

1from typing import Union 

2import numpy as np 

3import scipy 

4import scipy.constants 

5 

6 

7class HestonModel: 

8 def __init__(self, long_run_variance, mean_reversion_speed, vol_of_vol, 

9 initial_variance, correlation): 

10 """_summary_ 

11 

12 Args: 

13 long_run_variance (_type_): _description_ 

14 mean_reversion_speed (_type_): _description_ 

15 vol_of_vol (_type_): _description_ 

16 initial_variance (_type_): _description_ 

17 correlation (_type_): _description_ 

18 """ 

19 self._long_run_variance = long_run_variance 

20 self._mean_reversion_speed = mean_reversion_speed 

21 self._vol_of_vol = vol_of_vol 

22 self._initial_variance = initial_variance 

23 self._correlation = correlation 

24 

25 def feller_condition(self): 

26 """Return True if the model parameter fulfill the Feller condition 

27 ..: 

28 

29 Returns: 

30 bool: True->Feller condition is fullfilled 

31 """ 

32 return 2*self._mean_reversion_speed*self._long_run_variance>self._vol_of_vol > 0 

33 

34 def get_initial_value(self)->np.ndarray: 

35 """Return the initial value (x0, v0) 

36 

37 Returns: 

38 np.ndarray: Initial value. 

39 """ 

40 return np.array([1.0, self._initial_variance]) 

41 

42 def _characteristic_func(self, xi, s0, v0, tau): 

43 """Characteristic function needed internally to compute call prices with analytic formula. 

44 """ 

45 ixi = 1j * xi 

46 d = np.sqrt((self._mean_reversion_speed - ixi * self._correlation * self._vol_of_vol)**2 

47 + self._vol_of_vol**2 * (ixi + xi**2)) 

48 g = (self._mean_reversion_speed - ixi * self._correlation * self._vol_of_vol - d) / (self._mean_reversion_speed - ixi * self._correlation * self._vol_of_vol + d) 

49 ee = np.exp(-d * tau) 

50 C = self._mean_reversion_speed * self._long_run_variance / self._vol_of_vol**2 * ( 

51 (self._mean_reversion_speed - ixi * self._correlation * self._vol_of_vol - d) * tau - 2. * np.log((1 - g * ee) / (1 - g)) 

52 ) 

53 D = (self._mean_reversion_speed - ixi * self._correlation * self._vol_of_vol - d) / self._vol_of_vol**2 * ( 

54 (1 - ee) / (1 - g * ee) 

55 ) 

56 return np.exp(C + D*v0 + ixi * np.log(s0)) 

57 

58 def call_price(self, s0: float, v0: float, K: Union[np.ndarray, float], ttm: Union[np.ndarray, float])->Union[np.ndarray, float]: 

59 """Computes a call price for the Heston model via integration over characteristic function. 

60 

61 Args: 

62 s0 (float): current spot 

63 v0 (float): current variance 

64 K (float): strike 

65 ttm (float): time to maturity 

66 """ 

67 if isinstance(ttm, np.ndarray): 

68 result = np.empty((ttm.shape[0], K.shape[0], )) 

69 for i in range(ttm.shape[0]): 

70 #for j in range(K.shape[0]): 

71 #result[i,j] = self.call_price(s0,v0,K[j], tau[i]) 

72 result[i,:] = self.call_price(s0,v0,K, ttm[i]) 

73 return result 

74 

75 def integ_func(xi, s0, v0, K, tau, num): 

76 ixi = 1j * xi 

77 if num == 1: 

78 return (self._characteristic_func(xi - 1j, s0, v0, tau) / (ixi * self._characteristic_func(-1j, s0, v0, tau)) * np.exp(-ixi * np.log(K))).real 

79 else: 

80 return (self._characteristic_func(xi, s0, v0, tau) / (ixi) * np.exp(-ixi * np.log(K))).real 

81 

82 if ttm < 1e-3: 

83 res = (s0-K > 0) * (s0-K) 

84 else: 

85 "Simplified form, with only one integration. " 

86 h = lambda xi: s0 * integ_func(xi, s0, v0, K, ttm, 1) - K * integ_func(xi, s0, v0, K, ttm, 2) 

87 res = 0.5 * (s0 - K) + 1/scipy.constants.pi * scipy.integrate.quad_vec(h, 0, 500.)[0] #vorher 500 

88 return res 

89 

90 

91 def apply_mc_step(self, x: np.ndarray, t0: float, t1: float, rnd: np.ndarray, inplace: bool = True, slv: np.ndarray= None): 

92 """Apply a MC-Euler step for the Heston Model for n different paths. 

93 

94 Args: 

95 x (np.ndarray): 2-d array containing the start values for the spot and variance. The first column contains the spot, the second the variance values. 

96 t0 (float): The current time. 

97 t1 (float): The next timestep to be computed. 

98 rnd (np.ndarray): Two-dimensional array of shape (n_sims,2) containing the normal random numbers. Each row of the array is used to compute the correlated random numbers for the respective simulation. 

99 slv (np.ndarray): Stochastic local variance (for each path) to be multiplied with the heston variance. This is used by the StochasticVolatilityModel and can be ignored. 

100 """ 

101 if not inplace: 

102 x_ = x.copy() 

103 else: 

104 x_ = x 

105 rnd_V = np.sqrt(1.0-self._correlation**2)*rnd[:,1] + self._correlation*rnd[:,0] 

106 rnd_corr_S = rnd[:,0] 

107 S = x_[:,0] 

108 v = x_[:,1] 

109 dt = t1-t0 

110 sqrt_dt = np.sqrt(dt) 

111 if slv is None: 

112 slv=1.0 

113 S *= np.exp(- 0.5*v*slv*dt + np.sqrt(v*slv)*rnd_corr_S*sqrt_dt) 

114 v += self._mean_reversion_speed*(self._long_run_variance-v)*dt + self._vol_of_vol*np.sqrt(v)*rnd_V*sqrt_dt 

115 x_[:,1] = np.maximum(v,0) 

116 return x_