Modeling-Short-Rate/chi_square_model.cpp at master · devDivij/Modeling-Short-Rate · GitHub

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
#include <iostream>
#include <random>
#include <vector>
using namespace std;

double normal_random_variable(int n, mt19937 &gen)
{
    normal_distribution<double> dist(0.0, 1.0);
    double z = dist(gen);
    return z;
}
double chi_square_random_variable(double d, mt19937 &gen)
{
    chi_squared_distribution<double> dist(d);
    double x = dist(gen);
    return x;
}
double poisson_random_variable(double lambda, mt19937 &gen)
{
    poisson_distribution<int> dist(lambda);
    double N = dist(gen);
    return N;
}
vector<double> short_rate(int n, double r0, double T, double volatility, double alpha, double b)
{
    double d = (4 * b * alpha) / pow(volatility, 2);
    vector<double> r(n + 1);
    double dt = T / n;

    r[0] = r0;

    double c = pow(volatility, 2) * (1 - exp(-alpha * dt)) / (4 * alpha);

    random_device rd;
    mt19937 gen(rd());
    if (d > 1)
    {
        for (int i = 0; i < n; ++i)
        {
            double Z = normal_random_variable(n, gen);
            double X = chi_square_random_variable(d - 1.0, gen);
            double lambda = r[i] * exp(-alpha * dt) / c;
            r[i + 1] = c * (pow((Z + sqrt(lambda)), 2) + X);
        }
    }
    else if (d <= 1)
    {

        for (int i = 0; i < n; ++i)
        {
            double lambda = r[i] * exp(-alpha * dt) / c;
            double N = poisson_random_variable(lambda / 2, gen);
            double X = chi_square_random_variable(d + 2 * N, gen);
            r[i + 1] = c * X;
        }
    }

    return r;
}