1337 lines
48 KiB
MQL5
1337 lines
48 KiB
MQL5
//+------------------------------------------------------------------+
|
|
//| distribution.mqh |
|
|
//| Copyright 2025, MetaQuotes Ltd. |
|
|
//| https://www.mql5.com |
|
|
//+------------------------------------------------------------------+
|
|
#property copyright "Copyright 2025, MetaQuotes Ltd."
|
|
#property link "https://www.mql5.com"
|
|
#include"base.mqh"
|
|
#include<Math\Stat\T.mqh>
|
|
#include<Math\Stat\Gamma.mqh>
|
|
#include"..\Utility\igami.mqh"
|
|
#include"..\Utility\igam.mqh"
|
|
//+------------------------------------------------------------------+
|
|
//| Constraints container structure for optimization |
|
|
//| Holds linear inequality constraint matrices (Matrix * X >= Vector)|
|
|
//+------------------------------------------------------------------+
|
|
struct Constraints
|
|
{
|
|
matrix _one; // Constraint coefficient mapping matrix (A)
|
|
vector _two; // Boundary value threshold vector (b)
|
|
|
|
// Default Constructor
|
|
Constraints(void)
|
|
{
|
|
_one = matrix::Zeros(0,0);
|
|
_two = vector::Zeros(0);
|
|
}
|
|
// Destructor
|
|
~Constraints(void)
|
|
{
|
|
}
|
|
// Copy Constructor
|
|
Constraints(Constraints &other)
|
|
{
|
|
_one = other._one;
|
|
_two = other._two;
|
|
}
|
|
// Overloaded Assignment Operator
|
|
void operator=(Constraints &other)
|
|
{
|
|
_one = other._one;
|
|
_two = other._two;
|
|
}
|
|
};
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Abstract base class defining the standard interface for continuous |
|
|
//| probability distributions used in conditional volatility models |
|
|
//+------------------------------------------------------------------+
|
|
class CDistribution
|
|
{
|
|
protected:
|
|
ulong m_num_params; // Total estimated parameters in the model
|
|
vector m_parameters; // Vector array housing parameter values
|
|
vector m_vector; // Multi-use workspace vector
|
|
string m_name; // Identifiable distribution title literal string
|
|
uint m_seed; // Active seed initialization flag state
|
|
CHighQualityRandStateShell m_generator; // Alglib high-quality random engine wrapper shell
|
|
ENUM_DISTRIBUTION_MODEL m_dist_model; // Target model classification identifier
|
|
bool m_initialized; // Status flag verifying functional configuration
|
|
|
|
// Internal sanity validator checking that parameters fall within valid logical bounds
|
|
vector _check_constraints(vector ¶ms_)
|
|
{
|
|
matrix bounds_ = bounds(vector::Zeros(0));
|
|
ulong nparams = params_.Size();
|
|
if(nparams != bounds_.Rows())
|
|
{
|
|
Print(__FUNCTION__, " error ", "parameters must have ", bounds_.Rows(), " elements");
|
|
return vector::Zeros(0);
|
|
}
|
|
if(!bounds_.Rows())
|
|
return vector::Zeros(0);
|
|
|
|
// Validate lower/upper boundaries sequentially: Matrix columns 0 and 1
|
|
for(ulong i = 0; i < bounds_.Rows(); ++i)
|
|
{
|
|
if(!(bounds_[i,0] <= params_[i] && params_[i] <= bounds_[i,1]))
|
|
{
|
|
Print(__FUNCTION__, " error ", ". Bounds violated ");
|
|
return vector::Zeros(0);
|
|
}
|
|
}
|
|
return params_;
|
|
}
|
|
|
|
// Set seed states across standard RTL or high-quality native math environments
|
|
bool _initialize_generator(int rstate, bool use_alglib=true)
|
|
{
|
|
if(use_alglib)
|
|
{
|
|
m_seed = uint(MathAbs(rstate));
|
|
|
|
if(rstate)
|
|
CHighQualityRand::HQRndSeed((int)m_seed, (int)m_seed+1, m_generator.GetInnerObj());
|
|
else
|
|
CHighQualityRand::HQRndRandomize(m_generator.GetInnerObj());
|
|
return true;
|
|
}
|
|
else
|
|
{
|
|
m_seed = uint(MathAbs(rstate));
|
|
MathSrand(m_seed);
|
|
}
|
|
return true;
|
|
}
|
|
|
|
// Subclass hook targeting parameterized setup sequences
|
|
virtual bool _initialize(vector& distribution_params, int seed=0)
|
|
{
|
|
return false;
|
|
}
|
|
|
|
public:
|
|
CDistribution(void):m_name(NULL),
|
|
m_dist_model(WRONG_VALUE),
|
|
m_num_params(0),
|
|
m_initialized(false),
|
|
m_parameters(vector::Zeros(0))
|
|
{
|
|
}
|
|
|
|
CDistribution(vector& params, int seed=0)
|
|
{
|
|
initialize(params, seed);
|
|
}
|
|
|
|
CDistribution(CDistribution& other)
|
|
{
|
|
m_name = other.name();
|
|
m_dist_model = other.distribution_type();
|
|
m_num_params = other.numParams();
|
|
m_parameters = other.get_params();
|
|
m_seed = other.randomState();
|
|
m_initialized = other.is_initialized();
|
|
}
|
|
|
|
void operator=(CDistribution& other)
|
|
{
|
|
m_name = other.name();
|
|
m_dist_model = other.distribution_type();
|
|
m_num_params = other.numParams();
|
|
m_parameters = other.get_params();
|
|
m_seed = other.randomState();
|
|
m_initialized = other.is_initialized();
|
|
}
|
|
|
|
~CDistribution(void)
|
|
{
|
|
}
|
|
|
|
// Prepares operational parameters and safety tracks before model evaluation loops
|
|
virtual bool initialize(vector& distribution_params, int seed=0)
|
|
{
|
|
if(is_initialized())
|
|
{
|
|
m_initialized = false;
|
|
}
|
|
return _initialize(distribution_params, seed);
|
|
}
|
|
|
|
virtual ENUM_DISTRIBUTION_MODEL distribution_type(void) { return m_dist_model; }
|
|
virtual bool is_initialized(void) { return m_initialized; }
|
|
virtual CHighQualityRandStateShell* generator(void) { return GetPointer(m_generator); }
|
|
virtual uint randomState(void) { return m_seed; }
|
|
virtual ulong numParams(void) { return m_num_params; }
|
|
virtual string name(void) { return m_name; }
|
|
virtual vector get_params(void) { return m_parameters; }
|
|
void set_params(vector ¶meters) { m_parameters = parameters; }
|
|
|
|
// Returns simulation matrices filled with generated random numbers
|
|
virtual matrix simulator(ulong nsize, ulong ncols = 1) { return matrix::Zeros(0,0); }
|
|
|
|
// Random Number Generation (RNG) outputs formatted to vector or matrix shapes
|
|
virtual vector rng(ulong nsize)
|
|
{
|
|
matrix out = simulator(nsize);
|
|
return out.Col(0);
|
|
}
|
|
virtual matrix rng(ulong nrows, ulong ncols)
|
|
{
|
|
return simulator(nrows, ncols);
|
|
}
|
|
|
|
// Default stubs for bounds, likelihood evaluation and distribution moments
|
|
virtual Constraints constraints(void) { return Constraints(); }
|
|
virtual matrix bounds(vector& resids) { return matrix::Zeros(0,0); }
|
|
|
|
// Computes log-likelihood vectors over residuals vectors given their local variance sigmas
|
|
virtual vector loglikelihood(vector ¶meters, vector& resids, vector& sigmas, bool individial = false)
|
|
{
|
|
return vector::Zeros(0);
|
|
}
|
|
virtual vector startingValues(vector &resids) { return vector::Zeros(0); }
|
|
virtual vector ppf(vector &pits, vector ¶meters) { return vector::Zeros(0); } // Percent Point Function (Inverse CDF)
|
|
virtual vector cdf(vector &resids, vector ¶meters) { return vector::Zeros(0); } // Cumulative Distribution Function
|
|
virtual double moment(int n, vector& parameters) { return EMPTY_VALUE; }
|
|
virtual double partialMoment(int n, vector ¶meters, double z = 0.0) { return EMPTY_VALUE; } // Integrated partial density bounds
|
|
};
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Standard Normal Distribution implementation (Gaussian) |
|
|
//+------------------------------------------------------------------+
|
|
class CNormal: public CDistribution
|
|
{
|
|
private:
|
|
// Double factorial function evaluating product reductions: k!! = k * (k-2) * (k-4)...
|
|
int factorial2(int k)
|
|
{
|
|
int out = k;
|
|
for(int i = k-2; i > 1; out*=i, i-=2);
|
|
return out;
|
|
}
|
|
|
|
// Formulates non-central moments centered at zero for standardized normal terms
|
|
double munp(ulong n)
|
|
{
|
|
if(!n)
|
|
return 1.;
|
|
if(MathMod(double(n), 2.) == 0)
|
|
return (double)factorial2(int(n)-1); // Even orders match double factorials
|
|
else
|
|
return 0.; // Odd standard symmetric distributions track zero
|
|
}
|
|
|
|
// Maps structural moments directly from standard parameters (Mean, Variance, Skewness, Kurtosis)
|
|
vector moment_from_stats(ulong n, double mu, double mu2, double g1, double g2)
|
|
{
|
|
vector out(1);
|
|
if(n == 0) { out[0] = 1.0; return out; }
|
|
else if(n == 1) { out[0] = (!mu) ? munp(1) : mu; }
|
|
else if(n == 2) { out[0] = (!mu2 || !mu) ? munp(2) : mu2 + mu * mu; }
|
|
else if(n == 3)
|
|
{
|
|
if(!g1 || !mu2) out[0] = munp(3);
|
|
else
|
|
{
|
|
double mu3 = g1 * pow(mu2, 1.5);
|
|
out[1] = mu3 + 3.0 * mu * mu2 + mu * mu * mu;
|
|
}
|
|
}
|
|
else if(n == 4)
|
|
{
|
|
if(!g1 || !g2) out[0] = munp(4);
|
|
else
|
|
{
|
|
double mu3, mu4;
|
|
mu4 = (g2 + 3.0) * pow(mu2, 2.0);
|
|
mu3 = g1 * pow(mu2, 1.5);
|
|
out[1] = mu4 + 4 * mu * mu3 + 6 * mu * mu * mu2 + mu * mu * mu * mu;
|
|
}
|
|
}
|
|
else out[1] = munp(n);
|
|
|
|
return out;
|
|
}
|
|
|
|
// Converts standardized normal moment evaluations across given location/scale anchors
|
|
vector norm_moment(ulong order, double loc = 0.0, double scale = 1.)
|
|
{
|
|
ulong n = order;
|
|
double mu = 0.0, mu1 = 1.0, g1 = 0.0, g2 = 0.0;
|
|
vector val(1);
|
|
val = moment_from_stats(n, mu, mu1, g1, g2);
|
|
val[0] = pow(scale, order) * val[0];
|
|
return val;
|
|
}
|
|
|
|
virtual bool _initialize(vector& distribution_params, int seed=0) override
|
|
{
|
|
m_name = "Normal";
|
|
m_parameters = distribution_params;
|
|
m_num_params = 0; // Gaussian has zero shape parameters to estimate
|
|
m_initialized = _initialize_generator(seed, true);
|
|
return m_initialized;
|
|
}
|
|
|
|
public:
|
|
CNormal(void) { m_dist_model = DIST_NORMAL; }
|
|
CNormal(vector& params, int seed=0) { initialize(params, seed); }
|
|
~CNormal(void) {}
|
|
|
|
virtual Constraints constraints(void) override { Constraints out; return out; }
|
|
virtual matrix bounds(vector &resids) { return matrix::Zeros(0,0); }
|
|
|
|
// Normal distribution log-likelihood calculation: -0.5 * (ln(2*pi) + ln(sigma) + e^2 / sigma)
|
|
virtual vector loglikelihood(vector ¶meters, vector& resids, vector& sigmas, bool individual = false) override
|
|
{
|
|
vector lls = -0.5 * (log(2 * M_PI) + log(sigmas) + pow(resids, 2.0) / sigmas);
|
|
if(individual)
|
|
return lls;
|
|
else
|
|
{
|
|
vector out(1);
|
|
out[0] = lls.Sum();
|
|
return out;
|
|
}
|
|
}
|
|
|
|
virtual vector startingValues(vector& resids) override { return vector::Zeros(0); }
|
|
|
|
// Simulates multiple series vectors via Matrix tracking interfaces
|
|
virtual matrix simulator(ulong nsize, ulong ncols = 1) override
|
|
{
|
|
CMatrixDouble out;
|
|
out.Resize(nsize, ncols);
|
|
CHighQualityRand::HQRndNormalM(m_generator.GetInnerObj(), int(nsize), int(ncols), out);
|
|
return out.ToMatrix();
|
|
}
|
|
|
|
// Inverse Normal CDF calculations map probabilities directly onto standardized values
|
|
virtual vector ppf(vector &pits, vector ¶meters) override
|
|
{
|
|
_check_constraints(parameters);
|
|
vector out(pits.Size());
|
|
for(ulong i = 0; i < out.Size(); ++i)
|
|
out[i] = CNormalDistr::InvNormalCDF(pits[i]);
|
|
return out;
|
|
}
|
|
|
|
// Computes cumulative density vectors using traditional mathematical error integrations
|
|
virtual vector cdf(vector &resids, vector ¶meters) override
|
|
{
|
|
set_params(parameters);
|
|
_check_constraints(parameters);
|
|
vector out(resids.Size());
|
|
for(ulong i = 0; i < out.Size(); out[i] = CNormalDistr::NormalCDF(resids[i]), ++i);
|
|
return out;
|
|
}
|
|
|
|
virtual double moment(int n, vector& parameters)
|
|
{
|
|
set_params(parameters);
|
|
if(n < 0) return double("nan");
|
|
vector out = norm_moment(ulong(n));
|
|
return out[0];
|
|
}
|
|
|
|
// Evaluates conditional distributions inside partial integration bounds using recursive derivations
|
|
virtual double partialMoment(int n, vector ¶meters, double z = 0.0)
|
|
{
|
|
set_params(parameters);
|
|
if(n < 0) return double("nan");
|
|
else if(n == 0) return CNormalDistr::NormalCDF(z);
|
|
else if(n == 1) return -1. * CNormalDistr::NormalPDF(z);
|
|
else return -pow(z, double(n-1)) * CNormalDistr::NormalPDF(z) + double(n-1) * partialMoment(n-2, parameters, z);
|
|
}
|
|
};
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Standardized Student's t-Distribution for fat-tailed residuals |
|
|
//+------------------------------------------------------------------+
|
|
class CStudentsT: public CDistribution
|
|
{
|
|
private:
|
|
// Maps empirical parameters to generalized moment profiles specifically tailored to Student's t domains
|
|
double moment_from_stats(int n, double _mu, double _mu2, double _g1, double _g2, double df, double loc=0., double scale = 1.)
|
|
{
|
|
double out;
|
|
if(n == 0) return 1.;
|
|
else if(n == 1) { out = (_mu == EMPTY_VALUE) ? double("nan") : _mu; }
|
|
else if(n == 2) { out = (_mu2 == EMPTY_VALUE || _mu == EMPTY_VALUE) ? double("nan") : _mu2 + _mu * _mu; }
|
|
else if(n == 3)
|
|
{
|
|
if(_g1 == EMPTY_VALUE || _mu2 == EMPTY_VALUE || _mu == EMPTY_VALUE) out = double("nan");
|
|
else
|
|
{
|
|
double mu3 = _g1 * pow(_mu2, 1.5);
|
|
out = mu3 + 3. * _mu * _mu2 + _mu * _mu * _mu;
|
|
}
|
|
}
|
|
else if(n == 4)
|
|
{
|
|
if(_g1 == EMPTY_VALUE || _g2 == EMPTY_VALUE || _mu2 == EMPTY_VALUE || _mu == EMPTY_VALUE) out = double("nan");
|
|
else
|
|
{
|
|
double mu4 = (_g2 + 3.) * pow(_mu2, 2.);
|
|
double mu3 = _g1 * pow(_mu2, 1.5);
|
|
out = mu4 + 4. * _mu * mu3 + 6. * _mu * _mu * _mu2 + _mu * pow(_mu, 3);
|
|
}
|
|
}
|
|
else out = double("nan");
|
|
|
|
return out;
|
|
}
|
|
|
|
// Tracks internal moment values depending directly on active degrees of freedom (df) thresholds
|
|
vector _moment(int order, double df, double loc = 0.0, double scale = 1.)
|
|
{
|
|
vector out(0);
|
|
ulong n = ulong(order);
|
|
bool i0 = (scale > 0), i1 = (i0 && loc == 0.0), i2 = (i0 && loc != 0.0);
|
|
|
|
if(order < 0)
|
|
{
|
|
Print(__FUNCTION__, " Moment must be positive ");
|
|
return vector::Zeros(0);
|
|
}
|
|
double mu = EMPTY_VALUE, mu2 = EMPTY_VALUE, g1 = EMPTY_VALUE, g2 = EMPTY_VALUE;
|
|
if(0 < n && n < 5)
|
|
{
|
|
mu = df > 1. ? 0.0 : double("inf");
|
|
if(df > 1. && df <= 2.0) mu2 = double("inf");
|
|
else if(df > 2. && MathClassify(df) == FP_NORMAL) mu2 = df / (df - 2.);
|
|
else if(MathClassify(df) == FP_INFINITE) mu2 = 1.;
|
|
|
|
g1 = (df > 3.) ? 0.0 : double("nan");
|
|
if(df > 2. && df <= 4.) g2 = double("inf");
|
|
if(df > 4. && MathClassify(df) == FP_NORMAL) g2 = 6. / (df - 4.);
|
|
else if(MathClassify(df) == FP_INFINITE) g2 = 0.;
|
|
}
|
|
|
|
vector val(1);
|
|
val[0] = moment_from_stats(int(n), mu, mu2, g1, g2, df, loc, scale);
|
|
if(!i0) val[0] = double("nan");
|
|
if(i1 && loc == 0.0) val[0] = pow(scale, n) * val[0];
|
|
if(i2)
|
|
{
|
|
double res = 0., fac = scale / loc, valk;
|
|
for(int i = 0; i < int(n); ++i)
|
|
{
|
|
valk = moment_from_stats(i, mu, mu2, g1, g2, loc, scale);
|
|
res += double(np::comb(n, i, true)) * pow(fac, i) * valk;
|
|
}
|
|
res += pow(fac, n) * val[0];
|
|
res *= pow(loc, n);
|
|
val[0] = res;
|
|
}
|
|
return val;
|
|
}
|
|
|
|
// Integrates a localized snippet of the lower partial t-distribution curve
|
|
double _ord_t_partial_moment(int n, double z, double nu)
|
|
{
|
|
double mom = double("nan");
|
|
int ecode = 0;
|
|
if(n < 0 || double(n) > nu)
|
|
{
|
|
Print(__FUNCTION__, " n out of bounds ");
|
|
return mom;
|
|
}
|
|
else if(n == 0)
|
|
{
|
|
mom = MathCumulativeDistributionT(z, nu, ecode);
|
|
if(ecode) { Print(__FUNCTION__, " error with cdf ", ecode); return double("nan"); }
|
|
}
|
|
else if(n == 1)
|
|
{
|
|
double c = CGammaFunc::GammaFunc(0.5 + (n + 1.)) / (sqrt(nu * M_PI) * CGammaFunc::GammaFunc(0.5 * nu));
|
|
double e = 0.5 * (nu + 1.);
|
|
mom = (0.5 * (c * nu) / (1. - e)) * (pow(1. + pow(z, 2.) / nu, 1. - e));
|
|
}
|
|
else
|
|
{
|
|
ecode = 0;
|
|
double t1 = pow(z, n - 1.) * (nu + pow(z, 2.)) * MathProbabilityDensityT(z, nu, ecode);
|
|
if(ecode) { Print(__FUNCTION__, " error with cdf ", ecode); return double("nan"); }
|
|
double t2 = (n - 1) * nu * _ord_t_partial_moment(n - 2, z, nu);
|
|
mom = (1. / (double(n) - nu)) * (t1 - t2);
|
|
}
|
|
return mom;
|
|
}
|
|
|
|
virtual bool _initialize(vector& distribution_params, int seed=0) override
|
|
{
|
|
m_name = "Student's T";
|
|
m_num_params = 1; // 1 shape parameter: degrees of freedom
|
|
m_parameters = distribution_params;
|
|
m_initialized = _initialize_generator(seed, false);
|
|
return m_initialized;
|
|
}
|
|
|
|
public:
|
|
CStudentsT(void) { m_dist_model = DIST_STUDENT; }
|
|
CStudentsT(vector& params, int seed=0) { initialize(params, seed); }
|
|
~CStudentsT(void) {}
|
|
|
|
// Optimization bounds requiring df > 2.05 to maintain valid variance definitions
|
|
virtual Constraints constraints(void) override
|
|
{
|
|
Constraints out;
|
|
matrix x = {{1.0}, {-1.}};
|
|
vector y = {2.05, -500.};
|
|
out._one = x; out._two = y;
|
|
return out;
|
|
}
|
|
|
|
virtual matrix bounds(vector &resids) override
|
|
{
|
|
matrix out = {{2.05}, {500.}};
|
|
return out;
|
|
}
|
|
|
|
// Student's t Log-Likelihood accounting for fat-tailed kurtosis parameterizations
|
|
virtual vector loglikelihood(vector ¶meters, vector& resids, vector& sigmas, bool individual = false) override
|
|
{
|
|
set_params(parameters);
|
|
double nu = m_parameters[0];
|
|
double ll = log(CGammaFunc::GammaFunc((nu + 1.0) / 2.0)) - log(CGammaFunc::GammaFunc((nu) / 2.0)) - log(M_PI * (nu - 2.0)) / 2.;
|
|
vector lls = ll - (0.5 * log(sigmas));
|
|
lls -= ((nu + 1.) / 2.) * (log(1. + pow(resids, 2.0) / (sigmas * (nu - 2.))));
|
|
if(individual)
|
|
return lls;
|
|
else
|
|
{
|
|
vector out(1);
|
|
out[0] = lls.Sum();
|
|
return out;
|
|
}
|
|
}
|
|
|
|
// Estimates robust starting parameters using the sample kurtosis metrics
|
|
virtual vector startingValues(vector& resids) override
|
|
{
|
|
double k = np::kurtosis(resids, false, true);
|
|
double temp = k > 3.75 ? (4.0 * k - 6.0) / (k - 3.0) : 12;
|
|
double sv = MathMax(temp, 4.0);
|
|
vector out(1); out[0] = sv;
|
|
return out;
|
|
}
|
|
|
|
// Samples from the distribution and standardizes variance down to unity scale limits
|
|
virtual matrix simulator(ulong nsize, ulong ncols = 1) override
|
|
{
|
|
matrix out(nsize * ncols, 1);
|
|
vector temp; double vals[];
|
|
double std_dev = sqrt(m_parameters[0] / (m_parameters[0] - 2.));
|
|
if(!MathRandomT(m_parameters[0], int(nsize * ncols), vals) || !temp.Assign(vals) ||
|
|
!out.Col(temp, 0) || !out.Reshape(nsize, ncols))
|
|
{
|
|
Print(__FUNCTION__, " error sampling from student distribution ", GetLastError());
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
return out / (std_dev);
|
|
}
|
|
|
|
virtual vector ppf(vector &pits, vector ¶meters) override
|
|
{
|
|
set_params(parameters);
|
|
_check_constraints(parameters);
|
|
vector out(pits.Size());
|
|
int ecode = 0;
|
|
double var = m_parameters[0] / (m_parameters[0] - 2.0);
|
|
double scale = 1.0 / sqrt(var);
|
|
for(ulong i = 0; i < out.Size(); ++i)
|
|
{
|
|
out[i] = MathQuantileT(pits[i], m_parameters[0], ecode) * scale;
|
|
if(ecode) { Print(__FUNCTION__, " ppf() error ", ecode); return vector::Zeros(0); }
|
|
}
|
|
return out;
|
|
}
|
|
|
|
virtual vector cdf(vector &resids, vector ¶meters) override
|
|
{
|
|
set_params(parameters);
|
|
_check_constraints(parameters);
|
|
vector out(resids.Size());
|
|
int ecode = 0;
|
|
double var = m_parameters[0] / (m_parameters[0] - 2.0);
|
|
double scale = 1. / sqrt(var);
|
|
vector y = resids / scale;
|
|
for(ulong i = 0; i < out.Size(); ++i)
|
|
{
|
|
out[i] = (MathCumulativeDistributionT(y[i], m_parameters[0], ecode));
|
|
if(ecode) { Print(__FUNCTION__, " cdf() error ", ecode); return vector::Zeros(0); }
|
|
}
|
|
return out;
|
|
}
|
|
|
|
virtual double moment(int n, vector& parameters)
|
|
{
|
|
set_params(parameters);
|
|
if(n < 0) return double("nan");
|
|
_check_constraints(parameters);
|
|
double nu = parameters[0];
|
|
double var = nu / (nu - 2.);
|
|
double scale = 1.0 / sqrt(var);
|
|
vector out = _moment(n, nu, 0., scale);
|
|
return out[0];
|
|
}
|
|
|
|
virtual double partialMoment(int n, vector ¶meters, double z = 0.0)
|
|
{
|
|
set_params(parameters);
|
|
double nu = parameters[0];
|
|
double var = nu / (nu - 2.);
|
|
double scale = 1. / sqrt(var);
|
|
return pow(scale, n) * _ord_t_partial_moment(n, z / scale, nu);
|
|
}
|
|
};
|
|
|
|
//+---------------------------------------------------------------------+
|
|
//| Standardized Skewed Student's t-Distribution (Hansen 1994 setup) |
|
|
//| Handles both heavy-tail kurtosis and asymmetric financial profiling |
|
|
//+---------------------------------------------------------------------+
|
|
class CSkewStudent: public CDistribution
|
|
{
|
|
private:
|
|
// Internal distribution constants configuration (a, b, c) tracking skew adjustments
|
|
double _const_c(vector& parameters)
|
|
{
|
|
double eta = parameters[0];
|
|
return log(CGammaFunc::GammaFunc((eta + 1.) / 2.)) - log(CGammaFunc::GammaFunc(eta / 2.)) - log(M_PI * (eta - 2.)) / 2.;
|
|
}
|
|
double _const_b(vector& parameters)
|
|
{
|
|
double lam = parameters[1];
|
|
double a = _const_a(parameters);
|
|
return pow(1. + 3. * pow(lam, 2.) - pow(a, 2.), 0.5);
|
|
}
|
|
double _const_a(vector& parameters)
|
|
{
|
|
double eta = parameters[0];
|
|
double lam = parameters[1];
|
|
double c = _const_c(parameters);
|
|
return (4. * lam * exp(c) * (eta - 2.) / (eta - 1.));
|
|
}
|
|
|
|
double _ord_t_partial_moment(int n, double z, double nu)
|
|
{
|
|
double mom = double("nan"); int ecode = 0;
|
|
if(n < 0 || double(n) > nu)
|
|
{
|
|
Print(__FUNCTION__, " n out of bounds "); return mom;
|
|
}
|
|
else if(n == 0)
|
|
{
|
|
mom = MathCumulativeDistributionT(z, nu, ecode);
|
|
if(ecode) { Print(__FUNCTION__, " error with cdf ", ecode); return double("nan"); }
|
|
}
|
|
else if(n == 1)
|
|
{
|
|
double c = CGammaFunc::GammaFunc(0.5 + (n + 1.)) / (sqrt(nu * M_PI) * CGammaFunc::GammaFunc(0.5 * nu));
|
|
double e = 0.5 * (nu + 1.);
|
|
mom = (0.5 * (c * nu) / (1. - e)) * (pow(1. + pow(z, 2.) / nu, 1. - e));
|
|
}
|
|
else
|
|
{
|
|
ecode = 0;
|
|
double t1 = pow(z, n - 1.) * (nu + pow(z, 2.)) * MathProbabilityDensityT(z, nu, ecode);
|
|
if(ecode) { Print(__FUNCTION__, " error with cdf ", ecode); return double("nan"); }
|
|
double t2 = (n - 1) * nu * _ord_t_partial_moment(n - 2, z, nu);
|
|
mom = (1. / (double(n) - nu)) * (t1 - t2);
|
|
}
|
|
return mom;
|
|
}
|
|
|
|
virtual bool _initialize(vector& distribution_params, int seed=0) override
|
|
{
|
|
m_name = "Standardized Skew Honest Student's T";
|
|
m_num_params = 2; // 2 shape parameters: Tail thickness (eta) and Skewness (lambda)
|
|
m_parameters = distribution_params;
|
|
m_initialized = _initialize_generator(seed, true);
|
|
return m_initialized;
|
|
}
|
|
|
|
public:
|
|
CSkewStudent(void) { m_dist_model = DIST_SKEW_STUDENT; }
|
|
CSkewStudent(vector& params, int seed=0) { initialize(params, seed); }
|
|
~CSkewStudent(void) {}
|
|
|
|
// Multi-variable optimization parameters restricting asymmetry constraints to [-1, 1] range limits
|
|
virtual Constraints constraints(void) override
|
|
{
|
|
Constraints out;
|
|
matrix x = {{1, 0}, {-1, 0}, {0, 1}, {0, -1}};
|
|
vector y = {2.05, -300.0, -1, -1};
|
|
out._one = x; out._two = y;
|
|
return out;
|
|
}
|
|
virtual matrix bounds(vector &resids) override
|
|
{
|
|
matrix out = {{-1., 2.05}, {1., 300.}};
|
|
return out;
|
|
}
|
|
|
|
// Complex two-sided conditional likelihood calculations based on sign thresholds
|
|
virtual vector loglikelihood(vector ¶meters, vector& resids, vector& sigmas, bool individual = false) override
|
|
{
|
|
set_params(parameters);
|
|
double eta = parameters[0];
|
|
double lam = parameters[1];
|
|
double const_c = _const_c(parameters);
|
|
double const_a = _const_a(parameters);
|
|
double const_b = _const_b(parameters);
|
|
|
|
resids = resids / pow(sigmas, 0.5);
|
|
vector lls = log(const_b) + const_c - log(sigmas) / 2.;
|
|
if(MathAbs(lam) >= 1.) lam = np::sign(lam) * (1. - 1e-6);
|
|
|
|
vector signs(resids.Size());
|
|
for(ulong i = 0; i < signs.Size(); signs[i] = np::sign(resids[i] + const_a / const_b), ++i);
|
|
vector llf_resid = pow((const_b * resids + const_a) / (1. + signs * lam), 2.);
|
|
lls -= (eta + 1.) / 2. * log(1. + llf_resid / (eta - 2.));
|
|
|
|
if(individual) return lls;
|
|
else
|
|
{
|
|
vector out(1); out[0] = lls.Sum(); return out;
|
|
}
|
|
}
|
|
|
|
virtual vector startingValues(vector& resids) override
|
|
{
|
|
double k = np::kurtosis(resids, false, true);
|
|
double temp = k > 3.75 ? (4.0 * k - 6.0) / (k - 3.0) : 12;
|
|
double sv = (k > 3.75) ? MathMax((4. * k - 6.) / (k - 3.), 4.) : MathMax(12., 4.);
|
|
vector out(2); out[0] = sv; out[1] = 0.0; // Start with symmetrical assumption
|
|
return out;
|
|
}
|
|
|
|
// Map uniform arrays across functional inverse maps to run simulated sets
|
|
virtual matrix simulator(ulong nsize, ulong ncols = 1) override
|
|
{
|
|
matrix out(nsize * ncols, 1); vector uniforms(out.Rows());
|
|
CHighQualityRand::HQRndContinuous(m_generator.GetInnerObj(), int(out.Rows()), uniforms);
|
|
vector _ppf = ppf(uniforms, m_parameters);
|
|
out.Col(_ppf, 0); out.Reshape(nsize, ncols);
|
|
return out;
|
|
}
|
|
|
|
virtual vector ppf(vector &pits, vector ¶meters) override
|
|
{
|
|
set_params(parameters); _check_constraints(parameters);
|
|
vector out(pits.Size());
|
|
double eta = m_parameters[0], lam = m_parameters[1];
|
|
double a = _const_a(parameters), b = _const_b(parameters);
|
|
|
|
vector icdf = vector::Ones(pits.Size()) * -999.99;
|
|
int ecode = 0;
|
|
for(ulong i = 0; i < pits.Size(); ++i)
|
|
{
|
|
// Choose calculation fork based on regional probability split bounds
|
|
if(pits[i] < (1. - lam) / 2.)
|
|
icdf[i] = MathQuantileT((pits[i] / (1. - lam)), eta, ecode);
|
|
else
|
|
icdf[i] = MathQuantileT(((pits[i] - (1. - lam) / 2) / (1. + lam)), eta, ecode);
|
|
icdf[i] = icdf[i] * (1. + np::sign(pits[i] - (1. - lam) / 2.) * lam) * pow(1. - 2. / eta, 0.5 - a);
|
|
}
|
|
icdf = icdf / b;
|
|
return icdf;
|
|
}
|
|
|
|
virtual vector cdf(vector &resids, vector ¶meters) override
|
|
{
|
|
set_params(parameters); _check_constraints(parameters);
|
|
double eta = m_parameters[0], lam = m_parameters[1];
|
|
double a = _const_a(parameters), b = _const_b(parameters);
|
|
double var = eta / (eta - 2.);
|
|
|
|
vector y1 = (b * resids + a) / (1. - lam) * sqrt(var);
|
|
vector y2 = (b * resids + a) / (1. + lam) * sqrt(var);
|
|
vector p = vector::Zeros(resids.Size());
|
|
int ecode = 0;
|
|
for(ulong i = 0; i < p.Size(); ++i)
|
|
{
|
|
if(resids[i] < (-a / b))
|
|
p[i] = (1. - lam) * MathCumulativeDistributionT(y1[i], eta, ecode);
|
|
else
|
|
p[i] += ((1. - lam) / 2. + (1 + lam) * (MathCumulativeDistributionT(y2[i], eta, ecode) - 0.5));
|
|
}
|
|
return p;
|
|
}
|
|
|
|
// Blends left-tail and right-tail asymmetric components over combinatoric tracking fields
|
|
virtual double moment(int n, vector& parameters)
|
|
{
|
|
set_params(parameters);
|
|
if(n < 0) return double("nan");
|
|
_check_constraints(parameters);
|
|
|
|
double eta = m_parameters[0], lam = m_parameters[1];
|
|
double a = _const_a(parameters), b = _const_b(parameters);
|
|
|
|
double loc = -a / b;
|
|
double lscale = sqrt(1. - 2. / eta) * (1. - lam) / b;
|
|
double rscale = sqrt(1. - 2. / eta) * (1. + lam) / b;
|
|
double rpmom, lpmom, lhs, rhs, mom = 0.;
|
|
for(int k = 0; k < (n + 1); ++k)
|
|
{
|
|
rpmom = (0.5 * (CGammaFunc::GammaFunc(0.5 * (k + 1)) * CGammaFunc::GammaFunc(0.5 * (eta - k)) * pow(eta, (0.5 * k))) / (sqrt(M_PI) * CGammaFunc::GammaFunc(0.5 * eta)));
|
|
lpmom = pow(-1, k) * rpmom;
|
|
lhs = (1 - lam) * pow(lscale, k) * pow(loc, (n - k)) * lpmom;
|
|
rhs = (1 + lam) * pow(rscale, k) * pow(loc, (n - k)) * rpmom;
|
|
mom += np::comb(ulong(n), double(k)) * (lhs + rhs);
|
|
}
|
|
return mom;
|
|
}
|
|
|
|
virtual double partialMoment(int n, vector ¶meters, double z = 0.0)
|
|
{
|
|
set_params(parameters);
|
|
if(n < 0 || double(n) >= m_parameters[0]) return double("nan");
|
|
_check_constraints(parameters);
|
|
|
|
double eta = m_parameters[0], lam = m_parameters[1];
|
|
double a = _const_a(parameters), b = _const_b(parameters);
|
|
|
|
double loc = -a / b;
|
|
double lscale = sqrt(1. - 2. / eta) * (1. - lam) / b;
|
|
double rscale = sqrt(1. - 2. / eta) * (1. + lam) / b;
|
|
double lbound, lhs, rhs, mom = 0.;
|
|
for(int k = 0; k < (n + 1); ++k)
|
|
{
|
|
lbound = MathMin(z, loc);
|
|
lhs = ((1 - lam) * pow(loc, (n - k)) * pow(lscale, k) * _ord_t_partial_moment(k, (lbound - loc) / lscale, eta));
|
|
if(z > loc)
|
|
rhs = ((1 + lam) * pow(loc, (n - k)) * pow(rscale, k) * (_ord_t_partial_moment(k, (z - loc) / rscale, eta) - _ord_t_partial_moment(k, 0.0, eta)));
|
|
else
|
|
rhs = 0.0;
|
|
mom += np::comb(ulong(n), double(k)) * (lhs + rhs);
|
|
}
|
|
return mom;
|
|
}
|
|
};
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Generalized Error distribution for use with ARCH models |
|
|
//+------------------------------------------------------------------+
|
|
class CGeneralizedError: public CDistribution
|
|
{
|
|
private:
|
|
//--- Calculates the standardized partial moment of the Generalized Normal distribution
|
|
double _ord_gennorm_partial_moment(int n, int z, double beta)
|
|
{
|
|
if(n<0)
|
|
return double("nan");
|
|
|
|
// Math constants and scaling factor using Gamma Function
|
|
double w = 0.5*beta/CGammaFunc::GammaFunc(1./beta);
|
|
|
|
// Compute left-side integration boundary using the regularized complementary upper incomplete gamma function (igamc)
|
|
double lz = pow(MathAbs(MathMin(z,0)),beta);
|
|
double lterm = (w* (pow((-1), n))* (1 / beta)* CGammaFunc::GammaFunc((n + 1) / beta)* special::cephes::igamc((n + 1) / beta, lz));
|
|
|
|
// Compute right-side integration boundary
|
|
double rz = pow(MathMax(0,z),beta);
|
|
double rterm = w * (1 / beta) * CGammaFunc::GammaFunc((n + 1) / beta) * special::cephes::igamc((n + 1) / beta, rz);
|
|
|
|
return lterm+rterm;
|
|
}
|
|
|
|
//--- Computes statistical properties (mean, variance, skewness, excess kurtosis) using shape parameter beta
|
|
void stats(double beta, double &a, double &b, double &c, double &d)
|
|
{
|
|
// Use log-gamma variants to avoid numeric overflow
|
|
double c1 = log(CGammaFunc::GammaFunc(1./beta));
|
|
double c3 = log(CGammaFunc::GammaFunc(3./beta));
|
|
double c5 = log(CGammaFunc::GammaFunc(5./beta));
|
|
|
|
a = 0.; // Mean (Always zero due to symmetry)
|
|
b = exp(c3-c1); // Variance
|
|
c = 0.; // Skewness (Symmetric distribution)
|
|
d = exp(c5+c1-2.*c3)-3.; // Excess Kurtosis
|
|
return;
|
|
}
|
|
|
|
//--- Converts central moments to raw moments given location, scale, and statistical traits
|
|
double moment_from_stats(int n, double _mu,double _mu2, double _g1, double _g2, double df, double loc=0., double scale = 1.)
|
|
{
|
|
double out;
|
|
if(n==0)
|
|
return 1.;
|
|
else
|
|
if(n==1)
|
|
{
|
|
if(_mu == EMPTY_VALUE)
|
|
out = double("nan");
|
|
else
|
|
out = _mu;
|
|
}
|
|
else
|
|
if(n==2)
|
|
{
|
|
if(_mu2 == EMPTY_VALUE || _mu == EMPTY_VALUE)
|
|
out = double("nan");
|
|
else
|
|
out = _mu2 + _mu*_mu;
|
|
}
|
|
else
|
|
if(n==3)
|
|
{
|
|
if(_g1 == EMPTY_VALUE|| _mu2 == EMPTY_VALUE || _mu == EMPTY_VALUE)
|
|
out = double("nan");
|
|
else
|
|
{
|
|
double mu3 = _g1*pow(_mu2,1.5);
|
|
out = mu3+3.*_mu*_mu2+_mu*_mu*_mu;
|
|
}
|
|
}
|
|
else
|
|
if(n==4)
|
|
{
|
|
if(_g1 == EMPTY_VALUE||_g2 == EMPTY_VALUE||_mu2 == EMPTY_VALUE||_mu == EMPTY_VALUE)
|
|
out = double("nan");
|
|
else
|
|
{
|
|
double mu4 =(_g2+3.)*pow(_mu2,2.);
|
|
double mu3 = _g1*pow(_mu2,1.5);
|
|
out = mu4+4.*_mu*mu3+6.*_mu*_mu*_mu2+_mu*pow(_mu,3);
|
|
}
|
|
}
|
|
else
|
|
out = double("nan");
|
|
|
|
return out;
|
|
}
|
|
|
|
//--- Core implementation for raw nth-order central moment calculations
|
|
vector _moment(int order, double df,double loc = 0.0,double scale = 1.)
|
|
{
|
|
vector out(0);
|
|
ulong n = ulong(order);
|
|
bool i0,i1,i2;
|
|
i0 = i1 = i2 = false;
|
|
i0 = (scale>0);
|
|
i1 = (i0 & loc == 0.0);
|
|
i2 = (i0 & loc != 0.0);
|
|
|
|
if(order<0)
|
|
{
|
|
Print(__FUNCTION__, " Moment must be positive ");
|
|
return vector::Zeros(0);
|
|
}
|
|
|
|
double mu,mu2,g1,g2;
|
|
mu = mu2=g1 = g2 = EMPTY_VALUE;
|
|
|
|
if(0<n && n<5)
|
|
{
|
|
stats(df,mu,mu2,g1,g2);
|
|
}
|
|
|
|
vector val(1);
|
|
val[0] = moment_from_stats(int(n),mu,mu2,g1,g2,df,loc,scale);
|
|
|
|
if(!i0)
|
|
val[0] = double("nan");
|
|
if(i1 && loc == 0.0)
|
|
val[0] = pow(scale,n)*val[0];
|
|
|
|
// Apply binomial expansions if location shifting is required (loc != 0.0)
|
|
if(i2)
|
|
{
|
|
vector mom(4);
|
|
mom[0] = mu; mom[1] = mu2; mom[2] = g1; mom[3] = g2;
|
|
double res = 0.;
|
|
double fac = scale/loc;
|
|
double valk;
|
|
for(int i = 0; i<int(n); ++i)
|
|
{
|
|
valk = moment_from_stats(i,mu,mu2,g1,g2,loc,scale);
|
|
res += np::comb(n,i,true)*pow(fac,i)*valk;
|
|
}
|
|
res+=pow(fac,n)*val[0];
|
|
res*=pow(loc,n);
|
|
val[0] = res;
|
|
}
|
|
return val;
|
|
}
|
|
|
|
//--- Percent Point Function (Inverse CDF) using Cephes inverse incomplete gamma logic
|
|
double _ppf(double q, double beta)
|
|
{
|
|
double c = np::sign(q - 0.5);
|
|
return c*pow(special::cephes::igamci(1./beta,(1.0 + c) - 2.0*c*q),1./beta);
|
|
}
|
|
|
|
//--- Wrapper for scalar PPF calls with out-of-bounds safety checks
|
|
double gennorm_ppf(double nu,double loc, double scale, double q)
|
|
{
|
|
if(q<=0 || q>=1 || scale<=0)
|
|
{
|
|
double lb = -double("inf");
|
|
double ub = double("inf");
|
|
if(q == 0.) return lb;
|
|
else if(q==1.) return ub;
|
|
}
|
|
return _ppf(q,nu)*scale+loc;
|
|
}
|
|
|
|
//--- Vectorized version of the Percent Point Function
|
|
vector gennorm_ppfv(double nu,double loc, double scale, vector& q)
|
|
{
|
|
vector out(q.Size());
|
|
for(ulong i = 0; i<out.Size(); ++i)
|
|
out[i] = gennorm_ppf(nu,loc,scale,q[i]);
|
|
return out;
|
|
}
|
|
|
|
//--- Private override initialization sequence called via base class interface
|
|
virtual bool _initialize(vector& distribution_params, int seed=0) override
|
|
{
|
|
m_name = "Generalized Error Distribution";
|
|
m_num_params = 1;
|
|
m_parameters = distribution_params;
|
|
m_initialized= _initialize_generator(seed,false);
|
|
return m_initialized;
|
|
}
|
|
|
|
public:
|
|
//+------------------------------------------------------------------+
|
|
//| Constructors & Destructor |
|
|
//+------------------------------------------------------------------+
|
|
CGeneralizedError(void)
|
|
{
|
|
m_dist_model = DIST_GEN_ERROR;
|
|
}
|
|
|
|
CGeneralizedError(vector& params, int seed=0)
|
|
{
|
|
initialize(params,seed);
|
|
}
|
|
|
|
~CGeneralizedError(void)
|
|
{
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Parameter Constraints & Optimization Bounds |
|
|
//+------------------------------------------------------------------+
|
|
virtual Constraints constraints(void) override
|
|
{
|
|
Constraints out;
|
|
matrix x = {{1}, {-1}}; // Constraint definition matrix
|
|
vector y = {1.01, -500.0}; // Shape bounds (typically nu > 1.0)
|
|
out._one = x;
|
|
out._two = y;
|
|
return out;
|
|
}
|
|
|
|
virtual matrix bounds(vector &resids) override
|
|
{
|
|
matrix out = {{1.01},{500.}}; // Box constraints for optimizer [Min, Max]
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Log-Likelihood: Used during ARCH parameter estimation |
|
|
//+------------------------------------------------------------------+
|
|
virtual vector loglikelihood(vector ¶meters,vector& resids, vector& sigmas,bool individual = false) override
|
|
{
|
|
set_params(parameters);
|
|
double nu = parameters[0]; // Shape parameter (Degrees of freedom)
|
|
|
|
// Scaling constants for GED normalization
|
|
double logc = 0.5 * (-2. / nu * log(2.) + log(CGammaFunc::GammaFunc(1. / nu)) - log(CGammaFunc::GammaFunc(3. / nu)));
|
|
double c = exp(logc);
|
|
|
|
// Constant component of the log-likelihood function
|
|
double lls = log(nu) - logc - log(CGammaFunc::GammaFunc(1./nu)) - (1.+1./nu)*log(2.);
|
|
|
|
// Apply localized dynamic variance offsets over standard conditional variance (sigmas)
|
|
vector llsv = lls - (0.5 * log(sigmas));
|
|
llsv = llsv - (0.5*pow(MathAbs(resids/(sqrt(sigmas)*c)),nu));
|
|
|
|
if(individual)
|
|
return llsv; // Return pointwise vector for diagnostics
|
|
else
|
|
{
|
|
vector out(1);
|
|
out[0] = llsv.Sum(); // Return aggregated objective function value
|
|
return out;
|
|
}
|
|
}
|
|
|
|
//--- Optimization structural starting value guess
|
|
virtual vector startingValues(vector& resids) override
|
|
{
|
|
vector out(1);
|
|
out[0] = 1.5; // Standard starting shape value for financial residual optimization
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Simulator: Generates random matrix pathways via Gamma conversion |
|
|
//+------------------------------------------------------------------+
|
|
virtual matrix simulator(ulong nsize,ulong ncols = 1) override
|
|
{
|
|
double nu = m_parameters[0];
|
|
double gammas[];
|
|
double scale = sqrt(CGammaFunc::GammaFunc(3./nu)/CGammaFunc::GammaFunc(1./nu));
|
|
double shape = 1./nu;
|
|
|
|
// Generate standard gamma distributions to transform into GED random paths
|
|
if(!MathRandomGamma(shape,scale,int(nsize*ncols),gammas))
|
|
{
|
|
Print(__FUNCTION__," math gamma error ", GetLastError());
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
|
|
vector gammasv(nsize*ncols);
|
|
gammasv.Assign(gammas);
|
|
ArrayFree(gammas);
|
|
|
|
// Randomize direction signs uniformly to guarantee distribution symmetry
|
|
for(ulong i = 0; i<(nsize*ncols); gammasv[i]*= (2.*(double)CHighQualityRand::HQRndUniformI(m_generator.GetInnerObj(),2)-1.), ++i);
|
|
matrix out(nsize*ncols,1);
|
|
|
|
out.Col(gammasv,0);
|
|
out.Reshape(nsize,ncols);
|
|
return out/scale; // Rescale back to unit variance base
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Public Distribution Interfaces: PPF, CDF, Moments |
|
|
//+------------------------------------------------------------------+
|
|
virtual vector ppf(vector &pits, vector ¶meters) override
|
|
{
|
|
set_params(parameters);
|
|
_check_constraints(parameters);
|
|
double nu = m_parameters[0];
|
|
double c1 = log(CGammaFunc::GammaFunc(1./nu));
|
|
double c3 = log(CGammaFunc::GammaFunc(3./nu));
|
|
|
|
double var = exp(c3-c1);
|
|
double scale = 1./sqrt(var);
|
|
|
|
return gennorm_ppfv(nu,0.0,scale,pits);
|
|
}
|
|
|
|
virtual vector cdf(vector &resids, vector ¶meters) override
|
|
{
|
|
set_params(parameters);
|
|
_check_constraints(parameters);
|
|
|
|
double nu = m_parameters[0];
|
|
double c1 = log(CGammaFunc::GammaFunc(1./nu));
|
|
double c3 = log(CGammaFunc::GammaFunc(3./nu));
|
|
|
|
double var = exp(c3-c1);
|
|
double scale = 1./sqrt(var);
|
|
|
|
vector out = resids/scale;
|
|
|
|
// Classify outputs and safely compute localized probability values via igamc
|
|
for(ulong i = 0; i<out.Size(); ++i)
|
|
{
|
|
switch(MathClassify(out[i]))
|
|
{
|
|
case FP_INFINITE:
|
|
out[i] = 1.;
|
|
break;
|
|
case FP_NAN:
|
|
out[i] = double("nan");
|
|
break;
|
|
default:
|
|
out[i] = (0.5 + (0.5 * np::sign(out[i]))) - (0.5 * np::sign(out[i])) * special::cephes::igamc(1.0/nu, pow(MathAbs(out[i]),nu));
|
|
break;
|
|
}
|
|
}
|
|
return out;
|
|
}
|
|
|
|
virtual double moment(int n, vector& parameters)
|
|
{
|
|
set_params(parameters);
|
|
_check_constraints(parameters);
|
|
|
|
double nu = m_parameters[0];
|
|
double c1 = log(CGammaFunc::GammaFunc(1./nu));
|
|
double c3 = log(CGammaFunc::GammaFunc(3./nu));
|
|
|
|
double var = exp(c3-c1);
|
|
double scale = 1./sqrt(var);
|
|
|
|
vector out = _moment(n,nu,0.0,scale);
|
|
return out[0];
|
|
}
|
|
|
|
virtual double partialMoment(int n, vector ¶meters, double z = 0.0)
|
|
{
|
|
set_params(parameters);
|
|
_check_constraints(parameters);
|
|
|
|
double nu = m_parameters[0];
|
|
double c1 = log(CGammaFunc::GammaFunc(1./nu));
|
|
double c3 = log(CGammaFunc::GammaFunc(3./nu));
|
|
|
|
double var = exp(c3-c1);
|
|
double scale = 1./sqrt(var);
|
|
|
|
// Compute partial localized moment translated against current bounds scale
|
|
double mom = pow(scale,n)*_ord_gennorm_partial_moment(n,int(z/scale),nu);
|
|
return mom;
|
|
}
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| RNG |
|
|
//+------------------------------------------------------------------+
|
|
class BootstrapRng
|
|
{
|
|
protected:
|
|
vector m_std_resid; // Vector storing standardized residuals for empirical bootstrapping
|
|
ulong m_start; // Initial starting index for empirical sampling
|
|
ulong m_index; // Current index position tracked during random generation
|
|
CHighQualityRandStateShell m_rstate;// State wrapper for the High Quality Random Number Generator
|
|
CDistribution* m_distribution; // Pointer to the chosen parametric distribution model
|
|
bool m_initialized; // Flag tracking whether the generator is successfully ready for use
|
|
|
|
public:
|
|
//+------------------------------------------------------------------+
|
|
//| Default Constructor: Initializes pointers and states to safe defaults|
|
|
//+------------------------------------------------------------------+
|
|
BootstrapRng(void)
|
|
{
|
|
//--- m_rstate = NULL;
|
|
m_distribution = NULL;
|
|
m_initialized = false;
|
|
m_start = m_index = ULONG_MAX;
|
|
m_std_resid = vector::Zeros(0);
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Parametric Constructor: Instantiates an RNG using a specific |
|
|
//| mathematical distribution model. |
|
|
//+------------------------------------------------------------------+
|
|
BootstrapRng(ENUM_DISTRIBUTION_MODEL type, vector& parameters, int seed)
|
|
{
|
|
m_distribution = NULL;
|
|
|
|
// Instantiate the specific parametric distribution class based on type
|
|
switch(type)
|
|
{
|
|
case DIST_NORMAL:
|
|
m_distribution = new CNormal(parameters,seed);
|
|
break;
|
|
case DIST_STUDENT:
|
|
m_distribution = new CStudentsT(parameters,seed);
|
|
break;
|
|
case DIST_SKEW_STUDENT:
|
|
m_distribution = new CSkewStudent(parameters,seed);
|
|
break;
|
|
case DIST_GEN_ERROR:
|
|
m_distribution = new CGeneralizedError(parameters,seed);
|
|
break;
|
|
default:
|
|
Print(__FUNCTION__, " critical error. Failed instantiation of distribution ");
|
|
break;
|
|
}
|
|
|
|
// Verify successful initialization of the distribution object
|
|
m_initialized = (m_distribution.is_initialized()==true);//((CheckPointer(m_distribution)!=POINTER_INVALID) && (m_distribution.initialize(parameters,seed)));
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Empirical Constructor: Instantiates an RNG for historic |
|
|
//| bootstrapping using a fixed dataset of standardized residuals. |
|
|
//+------------------------------------------------------------------+
|
|
BootstrapRng(vector& std_resid, ulong start, int seed = 0)
|
|
{
|
|
m_initialized = false;
|
|
m_std_resid = std_resid;
|
|
m_start = start;
|
|
m_index = m_start;
|
|
m_distribution = NULL;
|
|
|
|
// Seed the High Quality Random Number Generator
|
|
if(seed)
|
|
CHighQualityRand::HQRndSeed(seed,seed+1,m_rstate.GetInnerObj());
|
|
else
|
|
CHighQualityRand::HQRndRandomize(m_rstate.GetInnerObj());
|
|
|
|
m_initialized = true;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Destructor: Clean up dynamically allocated distribution memory |
|
|
//+------------------------------------------------------------------+
|
|
~BootstrapRng(void)
|
|
{
|
|
if(CheckPointer(m_distribution) == POINTER_DYNAMIC)
|
|
delete m_distribution;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Check if the instance is ready to generate random data |
|
|
//+------------------------------------------------------------------+
|
|
bool is_initialized(void)
|
|
{
|
|
return m_initialized;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Generates a 1D vector of random samples via empirical bootstrap |
|
|
//+------------------------------------------------------------------+
|
|
vector rng(ulong size)
|
|
{
|
|
// Check if we are using empirical data (residuals vector populated)
|
|
if(m_std_resid.Size())
|
|
{
|
|
// Prevent out-of-bounds tracking if the pointer index exceeds available data
|
|
if(m_index>=m_std_resid.Size())
|
|
{
|
|
Print(__FUNCTION__, " not enough data points ");
|
|
return vector::Zeros(0);
|
|
}
|
|
|
|
// 1. Generate Uniform random values between [0, 1)
|
|
vector index = vector::Zeros(size);
|
|
for(ulong i = 0; i<size; index[i] = CHighQualityRand::HQRndUniformR(m_rstate.GetInnerObj()), ++i);
|
|
|
|
// 2. Map uniform samples to empirical indices within current bound constraints
|
|
index = floor(double(m_index)*index);
|
|
m_index+=1; // Increment current tracked horizon point
|
|
|
|
// 3. Extract the actual residual values from the historic vector matching the generated indices
|
|
vector out = index;
|
|
for(ulong i = 0; i<size; ++i)
|
|
out[i] = m_std_resid[ulong(index[i])];
|
|
|
|
return out;
|
|
}
|
|
return EMPTY_VECTOR;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Generates a 2D matrix of random samples using either parametric |
|
|
//| or empirical methods depending on how the class was instantiated.|
|
|
//+------------------------------------------------------------------+
|
|
matrix rng(ulong size, ulong cols)
|
|
{
|
|
// Route generation to parametric distribution if instantiated
|
|
if(CheckPointer(m_distribution)!=POINTER_INVALID)
|
|
return m_distribution.rng(size,cols);
|
|
else
|
|
{
|
|
// Otherwise, fallback to generating a 1D empirical vector and reshaping it into a matrix
|
|
matrix out(size*cols,1);
|
|
vector col = rng(size*cols);
|
|
out.Col(col,0);
|
|
out.Reshape(size,cols);
|
|
return out;
|
|
}
|
|
}
|
|
};
|
|
|