<commit_before>// <NAME> 2012.
//
// This example performs inference in the CoxIngersollRoss model
// (CIR model) which is a stochastic differential equation of the
// form
//
//    dr(t) = kappa*(theta - r(t))*dt + sigma*sqrt(r(t))*dW(r).
//
// The goal is to estimate kappa, theta and sigma from observed data.
// The CIR model is a model for interest rates.
//
// This will be done in two steps. First, an approximate log-likelihood
// function will be computed with a discretization of the PDE equal to
// the observed time interval. Second, using the first estimate as an
// initial estimate, a Monte-Carlo simulation will be used to compute
// a better likelihood function.
//
// This example illustrates that automatic differentiation can be used
// to compute first and second order derivatives of very complicated
// objective functions.
//
// The parameters optimized over are (log(kappa), log(theta), log(sigma)).
// This simple reparametrization prevents them from becoming negative.

#include <cmath>
#include <functional>
#include <iostream>
#include <random>

#include <spii/auto_diff_term.h>
#include <spii/solver.h>

using namespace spii;

// Probability density for a 1D normal distribution.
template<typename R>
R normpdf(const R& x, const R& mu, const R& sigma)
{
	R diff = (mu - x) / (sqrt(2.0)*sigma);
	return 1.0 / (sigma * sqrt(2.0 * 3.141592))
	       * exp(- diff * diff);
}

// Log probability density for a 1D normal distribution.
template<typename R>
R lognormpdf(const R& x, const R& mu, const R& sigma)
{
	R diff = (mu - x) / (sqrt(2.0)*sigma);
	return - log(sigma * sqrt(2.0 * 3.141592))
	       - diff * diff;
}

struct CIRApproximateLL
{
	const std::vector<double>& data;
	CIRApproximateLL(const std::vector<double>& data_in):
		data(data_in)
	{
	}

	template<typename R>
	R operator()(R* param) const
	{
		// Retrieve the three parameters.
		R kappa = exp(param[0]);
		R theta = exp(param[1]);
		R sigma  = exp(param[2]);

		R L = 0.0;
		for (size_t i = 0; i < data.size() - 1; ++i) {
			R eps_iplus1 = data[i+1] - data[i] - kappa*theta + kappa*data[i];
			L += lognormpdf<R>(eps_iplus1, 0, sigma * sqrt(data[i]));
		}

		// Return the negative log-likelihood since we are minimizing instead
		// of maximizing.
		return -L;
	}
};

struct SimulateCIR
{
	double r0;
	double r1;
	int M;
	int K;

	static unsigned seed;
	unsigned local_seed;

	SimulateCIR(double r0, double r1, int M, int K)
	{
		this->r0 = r0;
		this->r1 = r1;
		this->M = M;
		this->K = K;
		this->local_seed = seed++;
	}

	template<typename R>
	R operator()(R* param) const
	{
		// The random number generator
		std::seed_seq seed_s(&this->local_seed, &this->local_seed + 1);
		std::mt19937 prng(seed_s);
		std::normal_distribution<double> normal;
		auto randn = std::bind(normal, prng);

		// Retrieve the three parameters.
		R kappa = exp(param[0]);
		R theta = exp(param[1]);
		R sigma = exp(param[2]);

		double d = 1.0 / double(this->M);

		// Simulate K trajectories from r0 to r1
		R p = 0.0;
		for (int k = 1; k <= this->K; ++k) {
			R r = this->r0;
			// Simulate the trajectory in M-1 steps.
			for (int i = 1; i <= this->M - 1; ++i) {
				r = r + kappa*(theta - r)*d + sigma * sqrt(r*d) * randn();
				// This is to make sure r never becomes negative (which is
				// impossible in the continuous case).
				if (to_double(r) < 1e-8) {
					r = 1e-8;
				}
			}
			// Add the probability to go from r to r1.
			p += normpdf<R>(this->r1, r + kappa*(theta - r)*d, sqrt(d*r)*sigma);
		}
		// The estimated probability is the average of all trajectories.
		p /= double(K);

		// Return the negative log-likelihood since we are minimizing instead
		// of maximizing.
		return -log(p);
	}
};

unsigned SimulateCIR::seed = 0;

int main_function()
{
	// Read the observed data from standard input.
	std::vector<double> data;
	std::cerr << "Reading data from stdin... (use CIR.txt for example data)\n";
	while (true) {
		double r;
		std::cin >> r;
		if (! std::cin) {
			break;
		}
		data.push_back(r);
	}
	std::cerr << "Read " << data.size() << " data points\n";

	// These are the initial estimates of the parameters.
	// Using these with the simulated log-likelihood does
	// not work very well.
	double kappa = 1.0;
	double theta = 1.0;<commit_msg>Handle random number generation better.<commit_after>// <NAME> 2012.
//
// This example performs inference in the CoxIngersollRoss model
// (CIR model) which is a stochastic differential equation of the
// form
//
//    dr(t) = kappa*(theta - r(t))*dt + sigma*sqrt(r(t))*dW(r).
//
// The goal is to estimate kappa, theta and sigma from observed data.
// The CIR model is a model for interest rates.
//
// This will be done in two steps. First, an approximate log-likelihood
// function will be computed with a discretization of the PDE equal to
// the observed time interval. Second, using the first estimate as an
// initial estimate, a Monte-Carlo simulation will be used to compute
// a better likelihood function.
//
// This example illustrates that automatic differentiation can be used
// to compute first and second order derivatives of very complicated
// objective functions.
//
// The parameters optimized over are (log(kappa), log(theta), log(sigma)).
// This simple reparametrization prevents them from becoming negative.

#include <cmath>
#include <functional>
#include <iostream>
#include <random>

#include <spii/auto_diff_term.h>
#include <spii/solver.h>

using namespace spii;

// Probability density for a 1D normal distribution.
template<typename R>
R normpdf(const R& x, const R& mu, const R& sigma)
{
	R diff = (mu - x) / (sqrt(2.0)*sigma);
	return 1.0 / (sigma * sqrt(2.0 * 3.141592))
	       * exp(- diff * diff);
}

// Log probability density for a 1D normal distribution.
template<typename R>
R lognormpdf(const R& x, const R& mu, const R& sigma)
{
	R diff = (mu - x) / (sqrt(2.0)*sigma);
	return - log(sigma * sqrt(2.0 * 3.141592))
	       - diff * diff;
}

struct CIRApproximateLL
{
	const std::vector<double>& data;
	CIRApproximateLL(const std::vector<double>& data_in):
		data(data_in)
	{
	}

	template<typename R>
	R operator()(R* param) const
	{
		// Retrieve the three parameters.
		R kappa = exp(param[0]);
		R theta = exp(param[1]);
		R sigma  = exp(param[2]);

		R L = 0.0;
		for (size_t i = 0; i < data.size() - 1; ++i) {
			R eps_iplus1 = data[i+1] - data[i] - kappa*theta + kappa*data[i];
			L += lognormpdf<R>(eps_iplus1, 0, sigma * sqrt(data[i]));
		}

		// Return the negative log-likelihood since we are minimizing instead
		// of maximizing.
		return -L;
	}
};

typedef std::mt19937_64 Rng;

struct SimulateCIR
{
	double r0;
	double r1;
	int M;
	int K;

	static Rng::result_type seed;
	Rng::result_type local_seed;

	SimulateCIR(double r0, double r1, int M, int K)
	{
		this->r0 = r0;
		this->r1 = r1;
		this->M = M;
		this->K = K;
		this->local_seed = seed++;
	}

	template<typename R>
	R operator()(R* param) const
	{
		// The random number generator
		Rng prng(this->local_seed);
		std::normal_distribution<double> normal;
		auto randn = std::bind(normal, prng);

		// Retrieve the three parameters.
		R kappa = exp(param[0]);
		R theta = exp(param[1]);
		R sigma = exp(param[2]);

		double d = 1.0 / double(this->M);

		// Simulate K trajectories from r0 to r1
		R p = 0.0;
		for (int k = 1; k <= this->K; ++k) {
			R r = this->r0;
			// Simulate the trajectory in M-1 steps.
			for (int i = 1; i <= this->M - 1; ++i) {
				r = r + kappa*(theta - r)*d + sigma * sqrt(r*d) * randn();
				// This is to make sure r never becomes negative (which is
				// impossible in the continuous case).
				if (to_double(r) < 1e-8) {
					r = 1e-8;
				}
			}
			// Add the probability to go from r to r1.
			p += normpdf<R>(this->r1, r + kappa*(theta - r)*d, sqrt(d*r)*sigma);
		}
		// The estimated probability is the average of all trajectories.
		p /= double(K);

		// Return the negative log-likelihood since we are minimizing instead
		// of maximizing.
		return -log(p);
	}
};

Rng::result_type SimulateCIR::seed = Rng::default_seed;

int main_function()
{
	// Read the observed data from standard input.
	std::vector<double> data;
	std::cerr << "Reading data from stdin... (use CIR.txt for example data)\n";
	while (true) {
		double r;
		std::cin >> r;
		if (! std::cin) {
			break;
		}
		data.push_back(r);
	}
	std::cerr << "Read " << data.size() << " data points\n";

	// These are the initial estimates of the parameters.
	// Using these with the simulated log-likelihood does
	// not work very well.
	double kappa = 1.0;
	double theta = 1.0;
{'fraction_non_alphanumeric': 0.14319166373960895, 'fraction_numerical': 0.018733169418100926, 'mean_word_length': 3.134559535333979, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 64, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6571408', 'n_tokens_mistral': 3309, 'n_tokens_neox': 2971, 'n_words': 1022}