1369 lines
41 KiB
MQL5
1369 lines
41 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 |
|
|
//+------------------------------------------------------------------+
|
|
struct Constraints
|
|
{
|
|
matrix _one;
|
|
vector _two;
|
|
|
|
Constraints(void)
|
|
{
|
|
_one = matrix::Zeros(0,0);
|
|
_two = vector::Zeros(0);
|
|
}
|
|
~Constraints(void)
|
|
{
|
|
}
|
|
Constraints(Constraints &other)
|
|
{
|
|
_one = other._one;
|
|
_two = other._two;
|
|
}
|
|
void operator=(Constraints &other)
|
|
{
|
|
_one = other._one;
|
|
_two = other._two;
|
|
}
|
|
};
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| base class for distributions |
|
|
//+------------------------------------------------------------------+
|
|
class CDistribution
|
|
{
|
|
protected:
|
|
ulong m_num_params;
|
|
vector m_parameters;
|
|
vector m_vector;
|
|
string m_name;
|
|
uint m_seed;
|
|
CHighQualityRandStateShell m_generator;
|
|
ENUM_DISTRIBUTION_MODEL m_dist_model;
|
|
bool m_initialized;
|
|
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);
|
|
for(ulong i = 0; i<bounds_.Rows(); ++i)
|
|
{
|
|
if(!(bounds_[i,0]<=params_[i]<=bounds_[i,1]))
|
|
{
|
|
Print(__FUNCTION__, " error ", ". Bounds violated ");
|
|
return vector::Zeros(0);
|
|
}
|
|
}
|
|
return params_;
|
|
}
|
|
|
|
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;
|
|
}
|
|
|
|
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)
|
|
{
|
|
|
|
}
|
|
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;
|
|
}
|
|
virtual matrix simulator(ulong nsize,ulong ncols = 1)
|
|
{
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
virtual vector rng(ulong nsize)
|
|
{
|
|
matrix out = simulator(nsize);
|
|
return out.Col(0);
|
|
}
|
|
virtual matrix rng(ulong nrows,ulong ncols)
|
|
{
|
|
return simulator(nrows,ncols);
|
|
}
|
|
virtual Constraints constraints(void)
|
|
{
|
|
return Constraints();
|
|
}
|
|
virtual matrix bounds(vector& resids)
|
|
{
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
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);
|
|
}
|
|
virtual vector cdf(vector &resids, vector ¶meters)
|
|
{
|
|
return vector::Zeros(0);
|
|
}
|
|
virtual double moment(int n, vector& parameters)
|
|
{
|
|
return EMPTY_VALUE;
|
|
}
|
|
virtual double partialMoment(int n, vector ¶meters, double z = 0.0)
|
|
{
|
|
return EMPTY_VALUE;
|
|
}
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//|Standard normal distribution for use with ARCH models |
|
|
//+------------------------------------------------------------------+
|
|
class CNormal: public CDistribution
|
|
{
|
|
private:
|
|
|
|
int factorial2(int k)
|
|
{
|
|
int out = k;
|
|
for(int i = k-2; i>1; out*=i, i-=2);
|
|
return out;
|
|
}
|
|
|
|
double munp(ulong n)
|
|
{
|
|
if(!n)
|
|
return 1.;
|
|
if(MathMod(double(n),2.)==0)
|
|
return (double)factorial2(int(n)-1);
|
|
else
|
|
return 0.;
|
|
}
|
|
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)
|
|
{
|
|
if(!mu)
|
|
out[0] = munp(1);
|
|
else
|
|
out[0] = mu;
|
|
}
|
|
else
|
|
if(n==2)
|
|
{
|
|
if(!mu2 || !mu)
|
|
out[0] = munp(2);
|
|
else
|
|
out[0] = 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;
|
|
|
|
}
|
|
vector norm_moment(ulong order, double loc = 0.0, double scale = 1.)
|
|
{
|
|
ulong n = order;
|
|
bool i0 = (scale>0.);
|
|
bool i1 = bool(!loc);
|
|
bool i2 = bool(loc);
|
|
|
|
double mu,mu1,g1, g2;
|
|
mu = mu1 = g1 = g2 = 0.0;
|
|
if(n>0 && n<5)
|
|
{
|
|
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;
|
|
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);
|
|
}
|
|
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);
|
|
}
|
|
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();
|
|
}
|
|
virtual vector ppf(vector &pits, vector ¶meters) override
|
|
{
|
|
_check_constraints(parameters);
|
|
vector out(pits.Size());
|
|
for(ulong i = 0; i<out.Size(); CNormalDistr::InvNormalCDF(pits[i]), ++i);
|
|
return out;
|
|
}
|
|
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];
|
|
}
|
|
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 distribution for use with ARCH models |
|
|
//+------------------------------------------------------------------+
|
|
class CStudentsT: public CDistribution
|
|
{
|
|
private:
|
|
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;
|
|
}
|
|
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)
|
|
{
|
|
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)
|
|
{
|
|
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 += 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;
|
|
}
|
|
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;
|
|
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)
|
|
{
|
|
}
|
|
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;
|
|
}
|
|
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;
|
|
}
|
|
}
|
|
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;
|
|
}
|
|
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 distribution for use with ARCH models |
|
|
//+---------------------------------------------------------------------+
|
|
class CSkewStudent: public CDistribution
|
|
{
|
|
private:
|
|
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 Student's T";
|
|
m_num_params = 2;
|
|
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)
|
|
{
|
|
|
|
}
|
|
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;
|
|
}
|
|
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;
|
|
return out;
|
|
}
|
|
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];
|
|
double lam = m_parameters[1];
|
|
|
|
//double c = _const_c(parameters);
|
|
double a = _const_a(parameters);
|
|
double b = _const_b(parameters);
|
|
|
|
vector icdf = vector::Ones(pits.Size()) * -999.99;
|
|
vector signs = vector::Zeros(pits.Size());
|
|
int ecode = 0;
|
|
for(ulong i = 0; i<pits.Size(); ++i)
|
|
{
|
|
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];
|
|
double lam = m_parameters[1];
|
|
|
|
double a = _const_a(parameters);
|
|
double 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;
|
|
}
|
|
virtual double moment(int n, vector& parameters)
|
|
{
|
|
set_params(parameters);
|
|
if(n<0)
|
|
return double("nan");
|
|
_check_constraints(parameters);
|
|
|
|
double eta = m_parameters[0];
|
|
double lam = m_parameters[1];
|
|
|
|
double a = _const_a(parameters);
|
|
double 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];
|
|
double lam = m_parameters[1];
|
|
|
|
double a = _const_a(parameters);
|
|
double 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:
|
|
double _ord_gennorm_partial_moment(int n, int z, double beta)
|
|
{
|
|
if(n<0)
|
|
return double("nan");
|
|
|
|
double w = 0.5*beta/CGammaFunc::GammaFunc(1./beta);
|
|
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));
|
|
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;
|
|
}
|
|
void stats(double beta, double &a, double &b, double &c, double &d)
|
|
{
|
|
double c1 = log(CGammaFunc::GammaFunc(1./beta));
|
|
double c3 = log(CGammaFunc::GammaFunc(3./beta));
|
|
double c5 = log(CGammaFunc::GammaFunc(5./beta));
|
|
|
|
a = 0.;
|
|
b = exp(c3-c1);
|
|
c = 0.;
|
|
d = exp(c5+c1-2.*c3)-3.;
|
|
return;
|
|
}
|
|
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;
|
|
}
|
|
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];
|
|
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;
|
|
}
|
|
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);
|
|
}
|
|
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;
|
|
}
|
|
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;
|
|
}
|
|
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:
|
|
CGeneralizedError(void)
|
|
{
|
|
m_dist_model = DIST_GEN_ERROR;
|
|
}
|
|
CGeneralizedError(vector& params, int seed=0)
|
|
{
|
|
initialize(params,seed);
|
|
}
|
|
~CGeneralizedError(void)
|
|
{
|
|
|
|
}
|
|
virtual Constraints constraints(void) override
|
|
{
|
|
Constraints out;
|
|
matrix x = {{1}, {-1}};
|
|
vector y = {1.01, -500.0};
|
|
out._one = x;
|
|
out._two = y;
|
|
return out;
|
|
}
|
|
virtual matrix bounds(vector &resids) override
|
|
{
|
|
matrix out = {{1.01},{500.}};
|
|
return out;
|
|
}
|
|
virtual vector loglikelihood(vector ¶meters,vector& resids, vector& sigmas,bool individual = false) override
|
|
{
|
|
set_params(parameters);
|
|
double nu = parameters[0];
|
|
double logc = 0.5 * (-2. / nu * log(2.) + log(CGammaFunc::GammaFunc(1. / nu)) - log(CGammaFunc::GammaFunc(3. / nu)));
|
|
double c = exp(logc);
|
|
double lls = log(nu) - logc - log(CGammaFunc::GammaFunc(1./nu)) - (1.+1./nu)*log(2.);
|
|
vector llsv = lls - (0.5 * log(sigmas));
|
|
llsv = llsv - (0.5*pow(MathAbs(resids/(sqrt(sigmas)*c)),nu));
|
|
if(individual)
|
|
return llsv;
|
|
else
|
|
{
|
|
vector out(1);
|
|
out[0] = llsv.Sum();
|
|
return out;
|
|
}
|
|
}
|
|
virtual vector startingValues(vector& resids) override
|
|
{
|
|
vector out(1);
|
|
out[0] = 1.5;
|
|
return out;
|
|
}
|
|
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;
|
|
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);
|
|
|
|
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;//do i have to divide by scale
|
|
}
|
|
|
|
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;
|
|
|
|
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);
|
|
|
|
double mom = pow(scale,n)*_ord_gennorm_partial_moment(n,int(z/scale),nu);
|
|
return mom;
|
|
}
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//|RNG |
|
|
//+------------------------------------------------------------------+
|
|
class BootstrapRng
|
|
{
|
|
protected:
|
|
vector m_std_resid;
|
|
ulong m_start,m_index;
|
|
CHighQualityRandStateShell m_rstate;
|
|
CDistribution* m_distribution;
|
|
bool m_initialized;
|
|
public:
|
|
BootstrapRng(void)
|
|
{
|
|
//m_rstate = NULL;
|
|
m_distribution = NULL;
|
|
m_initialized = false;
|
|
m_start = m_index = ULONG_MAX;
|
|
m_std_resid = vector::Zeros(0);
|
|
}
|
|
BootstrapRng(ENUM_DISTRIBUTION_MODEL type, vector& parameters, int seed)
|
|
{
|
|
m_distribution = NULL;
|
|
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;
|
|
}
|
|
m_initialized = (m_distribution.is_initialized()==true);//((CheckPointer(m_distribution)!=POINTER_INVALID) && (m_distribution.initialize(parameters,seed)));
|
|
}
|
|
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;
|
|
|
|
if(seed)
|
|
CHighQualityRand::HQRndSeed(seed,seed+1,m_rstate.GetInnerObj());
|
|
else
|
|
CHighQualityRand::HQRndRandomize(m_rstate.GetInnerObj());
|
|
m_initialized = true;
|
|
}
|
|
~BootstrapRng(void)
|
|
{
|
|
if(CheckPointer(m_distribution) == POINTER_DYNAMIC)
|
|
delete m_distribution;
|
|
}
|
|
bool is_initialized(void)
|
|
{
|
|
return m_initialized;
|
|
}
|
|
vector rng(ulong size)
|
|
{
|
|
if(m_std_resid.Size())
|
|
{
|
|
if(m_index>=m_std_resid.Size())
|
|
{
|
|
Print(__FUNCTION__, " not enough data points ");
|
|
return vector::Zeros(0);
|
|
}
|
|
vector index = vector::Zeros(size);
|
|
for(ulong i = 0; i<size; index[i] = CHighQualityRand::HQRndUniformR(m_rstate.GetInnerObj()), ++i);
|
|
index = floor(double(m_index)*index);
|
|
m_index+=1;
|
|
vector out = index;
|
|
for(ulong i = 0; i<size; ++i)
|
|
out[i] = m_std_resid[ulong(index[i])];
|
|
return out;
|
|
}
|
|
return EMPTY_VECTOR;
|
|
}
|
|
matrix rng(ulong size, ulong cols)
|
|
{
|
|
if(CheckPointer(m_distribution)!=POINTER_INVALID)
|
|
return m_distribution.rng(size,cols);
|
|
else
|
|
{
|
|
matrix out(size*cols,1);
|
|
vector col = rng(size*cols);
|
|
out.Col(col,0);
|
|
out.Reshape(size,cols);
|
|
return out;
|
|
}
|
|
}
|
|
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
|
|
|