digifi 3.0.10

General purpose financial library and framework for financial modelling, portfolio optimisation, and asset pricing.
Documentation 
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

use ndarray::Array1;
use crate::error::DigiFiError;
use crate::statistics::ProbabilityDistribution;
use crate::statistics::continuous_distributions::NormalDistribution;
use crate::utilities::compare_len;
use crate::random_generators::{RandomGenerator, uniform_generators::FibonacciGenerator};


/// Implements the Accept-Reject Method, a Monte Carlo technique for generating random samples from a probability distribution.
/// 
/// # Input
/// - `f_x`: Target probability density function
/// - `g_x`: Proposal probability density function
/// - `y_sample`: Sample from the proposal distribution
/// - `m`: Constant such that m*g(x) >= f(x) for all x
/// - `uniform_sample`: Sample from a uniform distribution
/// 
/// # Output
/// - Sample array generated from the target distribution
/// 
/// # Errors
/// - Returns an error if the arrays produced by `f_x` and `g_x` are of different lengths.
/// - Returns an error if arrays produced by `f_x` and `g_x` are different in length to y_sample.
/// 
/// # Links
/// - Wikipedia: <https://en.wikipedia.org/wiki/Rejection_sampling>
/// - Original Source: N/A
pub fn accept_reject(
    f_x: &impl ProbabilityDistribution, g_x : &impl ProbabilityDistribution, y_sample: &Array1<f64>, m: f64, uniform_sample: &Array1<f64>
) -> Result<Array1<f64>, DigiFiError> {
    let g: Array1<f64> = g_x.pdf_iter(y_sample.iter())?;
    let f: Array1<f64> = f_x.pdf_iter(y_sample.iter())?;
    compare_len(&g.iter(), &f.iter(), "g", "f")?;
    compare_len(&g.iter(), &y_sample.iter(), "g", "sample")?;
    let x_sample: Vec<f64> = (0..y_sample.len()).into_iter().fold(vec![], |mut x_sample, i| {
        if g[i].is_finite() && f[i].is_finite() && (uniform_sample[i] * m * g[i]) <= f[i] {
            x_sample.push(y_sample[i]);
        }
        x_sample
    } );
    Ok(Array1::from_vec(x_sample))
}


/// Uses the Inverse Transform Method to generate random samples from a specified probability distribution function.
/// 
/// # Input
/// - `dist`: Probability distribution whose inverse cumulative distribution function is used to generate a sample
/// - `sample_size`: Number of random samples to generate
/// 
/// # Output
/// - Sample array generated from the specified distribution
/// 
/// # Latex
/// - X = F^{-1}(U)
/// 
/// # Links
/// - Wikipedia: <https://en.wikipedia.org/wiki/Inverse_transform_sampling>
/// - Original SOurce: N/A
pub fn inverse_transform(f_x: &impl ProbabilityDistribution, sample_size: usize) -> Result<Array1<f64>, DigiFiError> {
    let u: Array1<f64> = FibonacciGenerator::new_shuffle(sample_size)?.generate()?;
    f_x.inverse_cdf_iter(u.iter())
}


/// The Box-Muller algorithm transforms uniformly distributed samples into samples distributed according to the standard normal distribution.
/// 
/// # Input
/// - `uniform_sample_1`: First array of uniform samples
/// - `uniform_sample_2`: Second array of uniform samples
/// 
/// # Output
/// - Two arrays, each containing normal distributed samples
/// 
/// # Latex
/// -  Z_{0} = \\sqrt{-2ln(U_{1})} \\cos(2\\pi U_{2})
/// -  Z_{1} = \\sqrt{-2ln(U_{1})} \\sin(2\\pi U_{2})
/// 
/// # Links
/// - Wikipedia: <https://en.wikipedia.org/wiki/Box-Muller_transform>
/// - Original Source: <https://doi.org/10.1214%2Faoms%2F1177706645>
pub fn box_muller(uniform_sample_1: &Array1<f64>, uniform_sample_2: &Array1<f64>) -> Result<(Array1<f64>, Array1<f64>), DigiFiError> {
    compare_len(&uniform_sample_1.iter(), &uniform_sample_2.iter(), "uniform_sample_1", "uniform_sample_2")?;
    let n: usize = uniform_sample_1.len();
    let (z_0, z_1) = uniform_sample_1.iter().zip(uniform_sample_2.iter())
        .fold((Vec::with_capacity(n), Vec::with_capacity(n)), |(mut z_0, mut z_1), (us_1, us_2)| {
            let mult_1: f64 = if (us_1 == &0.0) || (us_1 == &1.0) { 0.0 } else { (-2.0 * us_1.ln()).sqrt() };
            z_0.push(mult_1 * (2.0 * std::f64::consts::PI * us_2).cos());
            z_1.push(mult_1 * (2.0 * std::f64::consts::PI * us_2).sin());
            (z_0, z_1)
        } );
    Ok((Array1::from_vec(z_0), Array1::from_vec(z_1)))
}


/// The Marsaglia polar method for generating standard normal random variables from uniformly distributed random numbers.
/// 
/// # Output
/// - A pair of standard normal random variables
/// 
/// # Latex
/// - S = W^{2}_{1} + W^{2}_{2}
/// - If S < 1, then:
/// - Z_{0} = W_{1} \\sqrt{\\frac{-2ln(S)}{S}}
/// - Z_{1} = W_{2} \\sqrt{\\frac{-2ln(S)}{S}}
/// 
/// # Links
/// - Wikipedia: <https://en.wikipedia.org/wiki/Marsaglia_polar_method>
/// - Original Source: <https://doi.org/10.1137%2F1006063>
pub fn marsaglia() -> Result<(f64, f64), DigiFiError> {
    loop {
        let u: Array1<f64> = 2.0 * FibonacciGenerator::new_shuffle(2)?.generate()? - 1.0;
        let s: f64 = u[0].powi(2) + u[1].powi(2);
        if s < 1.0 {
            let t: f64 = (-2.0 * s.ln() / s).sqrt();
            return Ok((u[0] * t, u[1] * t));
        }
    }
}


/// The Ziggurat algorithm is a method for generating random samples from a normal distribution.
///
/// Note: This version of Ziggurat algorithm does not implement a fallback algorithm for the tail.
/// 
/// # Input
/// - `x_guess`: Initial guess values
/// - `sample_size`: Number of random samples to generate
/// - `max_iterations`: Maximum number of iterations for the algorithm
/// 
/// # Output
/// - Sample arrays generated from the normal distribution
/// 
/// # Links
/// - Wikipedia: <https://en.wikipedia.org/wiki/Ziggurat_algorithm>
/// - Original Source: <https://doi.org/10.1145/1464291.1464310>
pub fn ziggurat(x_guess: &Array1<f64>, sample_size: usize, max_iterations: usize) -> Result<Array1<f64>, DigiFiError> {
    let normal_dist: NormalDistribution = NormalDistribution::build(0.0, 1.0)?;
    let y: Array1<f64> = normal_dist.pdf_iter(x_guess.iter())?;
    let mut z: Array1<f64> = Array1::from_vec(vec![1.0; sample_size]);
    for j in 0..sample_size {
        let u_1: Array1<f64> = 2.0 * FibonacciGenerator::new_shuffle(max_iterations)?.generate()? - 1.0;
        let mut i: usize = 0;
        'inner: while i < max_iterations {
            // Generates a random index and ensures that it is in the range [0, x_guess.len()-1]
            let rand_index: usize = (((x_guess.len() - 2) as f64) * FibonacciGenerator::new_shuffle(1)?.generate()?[0]).ceil() as usize;
            let x: f64 = u_1[rand_index] * x_guess[rand_index];
            if x.abs() < x_guess[rand_index + 1] {
                z[j] = x;
                break 'inner;
            }
            let u: f64 = FibonacciGenerator::new_shuffle(1)?.generate()?[0];
            let y_: f64 = y[rand_index] + u * (y[rand_index + 1] - y[rand_index]);
            let point: f64 = normal_dist.pdf(x)?;
            if y_ < point {
                z[j] = x;
                break 'inner;
            }
            i += 1;
        }
    }
    Ok(z)
}