blob: aa580368a49054d7a84ce7a23767b65b5104092f (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
// Hurst Exponent
//
// hurst_exponent
//
// Hurst exponent via R/S analysis
//
// H = slope of linear regression of log(R/S) vs log(n)
//
// Inputs: 1
// Outputs: 1
// Constants: 1
#include <math.h>
#include "ts.h"
typedef struct {
uint_t k;
} ts_hurst_exponent_params_t;
ts_parameter_t const ts_hurst_exponent_p = { .n_inputs = 1, .n_outputs = 1, .n_constants = 1, .params_size = sizeof(ts_hurst_exponent_params_t) };
uint_t ts_hurst_exponent_offset( ts_constant_t const *c )
{
( void )c;
return ( 0 );
}
int ts_hurst_exponent(
fp_t * const *y, uint_t y_n,
fp_t const * const *x, uint_t x_n,
uint_t n, ts_constant_t const *c
)
{
TS_PROLOGUE( ts_hurst_exponent_p );
if ( !c ) return -1;
uint_t k = ((ts_hurst_exponent_params_t *)c)->k;
if (k == 0 || k > n) return -1;
for ( uint_t i = 0; i < n; ++i )
{
if (i + 1 < k) {
y[0][i] = 0;
continue;
}
uint_t data_len = i + 1;
uint_t num_lags = 0;
fp_t sum_log_s = 0;
fp_t sum_log_rs = 0;
fp_t sum_log_s2 = 0;
fp_t sum_log_s_log_rs = 0;
for ( uint_t s = 2; s <= k && s <= data_len / 2; ++s )
{
uint_t ns = data_len / s;
if (ns == 0) continue;
fp_t total_log_rs = 0;
uint_t valid_segs = 0;
for ( uint_t seg = 0; seg < ns; ++seg )
{
uint_t start = seg * s;
uint_t end = start + s;
if (end > data_len) break;
fp_t mean = 0;
for ( uint_t p = start; p < end; ++p )
mean += x[0][p];
mean /= (fp_t)s;
fp_t cumsum = 0;
fp_t max_dev = 0, min_dev = 0;
for ( uint_t p = start; p < end; ++p )
{
cumsum += x[0][p] - mean;
if (cumsum > max_dev) max_dev = cumsum;
if (cumsum < min_dev) min_dev = cumsum;
}
fp_t r = max_dev - min_dev;
fp_t std = 0;
for ( uint_t p = start; p < end; ++p )
{
fp_t d = x[0][p] - mean;
std += d * d;
}
std = sqrt(std / (fp_t)s);
if (std > 0 && r > 0) {
total_log_rs += log(r / std);
valid_segs++;
}
}
if (valid_segs > 0) {
fp_t log_s = log((fp_t)s);
fp_t mean_log_rs = total_log_rs / (fp_t)valid_segs;
sum_log_s += log_s;
sum_log_rs += mean_log_rs;
sum_log_s2 += log_s * log_s;
sum_log_s_log_rs += log_s * mean_log_rs;
num_lags++;
}
}
if (num_lags >= 2) {
fp_t denom = (fp_t)num_lags * sum_log_s2 - sum_log_s * sum_log_s;
if (fabs(denom) > 1e-300) {
fp_t h = ((fp_t)num_lags * sum_log_s_log_rs - sum_log_s * sum_log_rs) / denom;
y[0][i] = h;
} else {
y[0][i] = 0;
}
} else {
y[0][i] = 0;
}
}
return 0;
}
|