Coverage for src/scrilla/analysis/integration.py: 0%

1# This file is part of scrilla: https://github.com/chinchalinchin/scrilla. 

2 

3# scrilla is free software: you can redistribute it and/or modify 

4# it under the terms of the GNU General Public License version 3 

5# as published by the Free Software Foundation. 

6 

7# scrilla is distributed in the hope that it will be useful, 

8# but WITHOUT ANY WARRANTY; without even the implied warranty of 

9# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

10# GNU General Public License for more details. 

11 

12# You should have received a copy of the GNU General Public License 

13# along with scrilla. If not, see <https://www.gnu.org/licenses/> 

14# or <https://github.com/chinchalinchin/scrilla/blob/develop/main/LICENSE>. 

15 

16""" 

17A module of functions for integrating stochastic differential equations through numerial techniques like Monte Carlo simulation. 

18""" 

19from typing import Callable 

20from scipy.stats import norm 

21from scipy.integrate import quad 

22from numpy import isinf, inf 

23import random 

24 

25from scrilla import settings 

26from scrilla.util import errors 

27import scrilla.util.outputter as outputter 

28 

29logger = outputter.Logger("scrilla.analysis.integration", settings.LOG_LEVEL) 

30 

31 

32def generate_random_walk(periods: int): 

33 return [norm.ppf(random.uniform(0, 1)) for i in range(periods)] 

34 

35# Condition: E(integral of vol ^2 dt) < inf 

36 

37 

38def verify_volatility_condition(volatility_function: Callable) -> bool: 

39 r""" 

40 Ito calculus requires the class of volatility functions that scale a simple Weiner process (i.e., a process with independent, identically distributed increments with mean 0 and variance 1) satisfies the following 

41 

42 $$ \int_{0}^{\infty} {\sigma (t)}^2 \,dt < \infty $$ 

43 

44 In order words, since the mean of an Weiner process increment is 0 (and thus \\(Var(X) = E(X^2)\\), this inequality imposes a condition on the process described by scaling the Weiner process by the volatility function: it must have a probability distribution that actually exists. This function returns a `bool` that signals whether or not the given function satisifies this condition. 

45 

46 Parameters 

47 ---------- 

48 1. **volatility_function** : ``Callable`` 

49 A function that describes the volatility of a process as a function of time measured. Must accept scalar input. 

50 """ 

51 integral = quad(func=lambda x: volatility_function(x)**2, a=0, b=inf) 

52 return not isinf(x=integral) 

53 

54# Remember forward increments! 

55 

56 

57def ito_integral(mean_function: Callable, volatilty_function: Callable, upper_bound: float, iterations: int = 1000): 

58 r""" 

59 Approximates the expectation of an Ito integral from 0 to `upper_bound` for a Weiner process described by a mean and volatility that are functions of time, i.e. 

60 

61 $$ E\{d X(t)\} = E \{ \int_{0}^{\infty} {\mu (t)} \, dt + \int_{0}^{\infty} {\sigma(t)} \, dW(t) \} $$ 

62 

63 Parameters 

64 ---------- 

65 1. **mean_function** : ``Callable`` 

66 A function that describes the mean of a process as a function of time. Must accept scalar input. 

67 2. **volatility_function** : ``Callable`` 

68 A function that describes the volatility of a process as a function of time. Must accept scalar input. 

69 3. **upper_bound** : ``float`` 

70 Upper limit for the integral. 

71 4. **iterations** : ``int`` 

72 The number of iterations in the Reimann-Stieltjes sum. 

73 """ 

74 if verify_volatility_condition(volatility_function=volatilty_function): 

75 # NOTE: Compute Ito Integral using forward increments due to the fact 

76 # the value at the end of the interval can be conditioned on 

77 # the value at the beginning of the interval. In other words, Ito 

78 # Processes are Martingales, ie. E(X2 | X1) = X1. 

79 time_delta = upper_bound / settings.ITO_STEPS 

80 

81 averaged_integral = 0 

82 for j in range(iterations): 

83 random_walk = generate_random_walk(settings.ITO_STEPS) 

84 this_integral, current_contribution = 0, 0 

85 current_contribution = 0 

86 for i in range(settings.ITO_STEPS): 

87 current_time = time_delta * i 

88 current_mean = mean_function(current_time) 

89 current_vol = volatilty_function(current_time) 

90 current_contribution = current_mean * \ 

91 time_delta + current_vol * random_walk[i] 

92 this_integral += current_contribution 

93 averaged_integral += this_integral/j 

94 

95 return averaged_integral 

96 raise errors.UnboundedIntegral( 

97 'E(Integral(volatility_function^2 dt) from 0 to infinity < infinity)')