Statistics and regression

Sections: Statistics · Percent Rank · R-Squared · Spearman Rank Correlation · Linear Regression · Fractional Differencing

Statistics

DownsideSemiDeviation

func (s Series) DownsideSemiDeviation(n int) Series

DownsideSemiDeviation returns the downside semi-deviation of s over a rolling window of n bars, the root mean square deviation from the window’s mean of the values below that mean. For the n values ending at bar t, with mu their mean, it is sigma- = sqrt(S- / N-) (Poster’s equation 2), where S- is the sum of (x - mu)^2 over the values x < mu and N- is their number. The inequality is strict: a value equal to mu is on neither side, so N+ + N- <= n, where N+ counts the values above mu. The divisor is the side’s own count N-, not n, so sigma- is the root mean square deviation of the values below the mean. Poster’s indicator deviation bands replace Bollinger’s single standard deviation with one deviation per side of the mean, sigma- below it and sigma+ above it (see UpsideSemiDeviation).

The input is any series, typically an indicator: Poster applies the bands to a fast SMA of the close, to the RSI and to the stochastic.

The result is window-local, each bar from its own window, so it does not depend on where the history starts. It costs O(n) per bar (one pass to scale the window, two to centre it and one scan of the sides) and one scratch slice of length n per call, besides the result, a new series of len(s) values; s is left unchanged.

The mean is the method’s own, computed as in StdDevSample: the window is anchored at its first value, scaled by a power of two and averaged with a compensated sum. It agrees with Average(n) to rounding, not bit for bit, and the method does not call Average. A value counts as below the mean when its scaled centred deviation is negative and as above it when the deviation is positive, by strict comparisons with 0 rather than by math.Signbit, so a deviation of -0, like one of +0, is on neither side. Membership agrees with exact arithmetic on windows of integers, or of multiples of a common power of two, whose span, counted in units of that power and multiplied by n, is below 2^52: there the differences from the first value and their sum are exact, and the rounding of the mean does not move any value to another side. Ties with a decimal mean are fragile: a value equal to the decimal mean may count on either side or on neither, on data with few decimals such as an RSI to two decimals or prices to the tick, because binary floating point cannot hold the decimal mean and no mean of the stored values, whether this one, Average(n), a naive sum or their exact mean, finds such ties every time. That can move N- by one and sigma- by a step, a property of the count normalisation, not of this method.

In exact arithmetic, other conventions differ from this one by a factor that changes from bar to bar with N-, not by a constant: the Markowitz semi-deviation divides by n and is sigma- times sqrt(N-/n), and the convention of Kahn and Stefek divides by n-1, which gives sigma- times sqrt(N-/(n-1)). Sortino’s downside deviation about a fixed target tau also divides by n but needs no kernel; for a series x it is the composition

x.Sub(tau).Min(0).Squared().Sum(n).Div(float64(n)).Sqrt().SetN(n-1, math.NaN())

in which SetN masks the start-up bars, which Sum fills from fewer values. The two sides add up to the variance, to rounding: N+ sigma+^2 + N- sigma-^2 = n StdDev(n)^2 = (n-1) StdDevSample(n)^2, with a value at the mean adding to neither side. On a window symmetric about its mean with no value at the mean the two sides are equal, bit for bit in the region above where membership agrees with exact arithmetic, and both equal StdDev(n) to rounding, so bands built on the two semi-deviations about Average(n) are then Bollinger’s to rounding. Negation swaps the sides exactly: s.Mul(-1).UpsideSemiDeviation(n) equals s.DownsideSemiDeviation(n) bit for bit, for every input. At n = 2 both sides are |s[t] - s[t-1]|/2 for t >= 1 while the difference does not overflow, and StdDev(2) bit for bit while that difference lies between about 3e-154 and 1.9e154 in magnitude, outside which StdDev’s running moments underflow or overflow; a window whose difference overflows, such as Series{math.MaxFloat64, -math.MaxFloat64}, is scaled first here and gives math.MaxFloat64 where StdDev(2) gives +Inf.

Poster’s bands are four lines about Average(n) of the indicator; the two oversold lines lie below it at k1 and k2 times sigma-, with k2 > k1, and UpsideSemiDeviation gives the whole recipe, its parameters and its warm-up.

The warm-up reads 0, not NaN: bars 0 to n-2 are 0 and the first valid index is n-1. As 0 is also the value of a flat window, trim the warm-up with Window(-(n-1)) or mask it with SetN(n-1, math.NaN()). A window length n < 2 or n > len(s) gives len(s) zeros, as for StdDev.

A NaN or an infinity at bar p gives NaN at every bar whose window holds it, bars max(p, n-1) to min(p+n-1, len(s)-1), and later windows recover; warm-up bars stay 0 even where the input is NaN. Values that are not finite are tested before flatness, so a window of equal infinities is NaN.

A flat window, one whose values are all equal (-0 and +0 count as equal), gives +0, decided exactly rather than through a rounded mean. The rule that an empty side gives +0 is a safeguard: a window that is not flat has values on both sides in practice, as its computed mean lies strictly between its least and greatest values. A side with values has a positive sum of squares, so its result is positive or, on underflow, +0: every zero result is +0.

The result lies between 0 and the largest |x - mu| of the window, to rounding. A finite window gives a finite value unless the exact value reaches math.MaxFloat64 plus half a unit in its last place, where binary64 rounds to +Inf, or comes within rounding of that, when it is +Inf: with n = len(s), Series{math.MaxFloat64, -math.MaxFloat64, 1}, whose exact downside value is math.MaxFloat64 plus 1/3, gives math.MaxFloat64. With n = len(s), Series{math.MaxFloat64, -math.MaxFloat64} gives math.MaxFloat64 on both sides, as a window whose differences overflow is scaled first. Series{math.MaxFloat64, -math.MaxFloat64, math.MaxFloat64}, with n = len(s), gives +Inf here, where the exact value is 4/3 of math.MaxFloat64, while UpsideSemiDeviation gives 1.1984620899082105e308, 2/3 of math.MaxFloat64. At the other end, a window that is not flat gives 0 only when the exact value is about 2.5e-324 or less, and with n = len(s), Series{0, 5e-324} gives 0 on both sides.

Where its membership agrees with exact arithmetic, the result is within 4 units in the last place plus 2^-50 times the window’s range of the exact value, as the references check, and within the 4 units alone on continuous data about a level and on random walks, though not in general.

Sources: Poster, “Indicator Deviation Bands: Extending The Bollinger Bands Concept”, Stocks & Commodities V42:12 (November 2024), pp. 14-17, for the definition; Kahn and Stefek, “Heat, Light, and Downside Risk”, chapter 27 of Grinold and Kahn, Advances in Active Portfolio Management (2019), for the n-1 convention.

Kurtosis

func (s Series) Kurtosis(n int) Series

Kurtosis returns the rolling population excess kurtosis of s over windows of n bars: g2 = m4 / m2^2 - 3 for the n values ending at bar t, where m_k = (1/n) sum (x - mean)^k is their k-th central moment about their own mean (a normal distribution has excess kurtosis 0). This is Natenberg’s normalised form, scipy’s kurtosis(bias=True) and R’s moments::kurtosis minus 3; Pearson’s kurtosis m4 / m2^2 (R’s moments::kurtosis, scipy’s fisher=False) is Kurtosis + 3. KurtosisSample is the sample form of Excel and pandas; on a window that is not flat, the conversion between them is KurtosisSample = ((n+1) Kurtosis + 6)(n-1)/((n-2)(n-3)). Zero crossings depend on the form: in exact arithmetic, KurtosisSample < 0 exactly when Kurtosis < -6/(n+1). Kaufman’s skewness and kurtosis formulas in Stocks & Commodities (September 2026) are non-standard rescalings.

Kurtosis is unitless. Apply it to returns or another stationary series: on price levels it mostly measures the shape of the trend. The sampling error is large on short windows; for normal data the standard error is about sqrt(24/n). To rounding, Kurtosis >= -2 and Kurtosis + 3 >= Skewness^2 + 1. Every window that is not flat gives exactly -2 at n = 2 and -1.5 (to rounding) at n = 3, so use n >= 4.

The warm-up reads 0, not NaN: bars 0 to n-2 are 0 and the first valid index is n-1. Since 0 is also the excess kurtosis of a normal distribution, trim (Window(-(n-1))) or mask (SetN(n-1, math.NaN())) the warm-up before using the values as data. The minimum n is 2: an n below 2 or above len(s) gives len(s) zeros, aligned with s. A flat window, whose values are all equal (-0 and +0 count as equal), gives +0, a sentinel rather than a measurement, as in Standardise and Correlation. A NaN or an infinity at bar p gives NaN at bars max(p, n-1) to min(p+n-1, len(s)-1), the windows that hold it, and later windows recover. Finite windows of any magnitude, up to math.MaxFloat64, give finite values.

Each bar is computed from its own window in O(n) time, so the value at a bar does not depend on where a slice of history starts.

Sources: Natenberg, Option Volatility and Pricing (2nd ed., 2015), Appendix B; Joanes and Gill, “Comparing measures of sample skewness and kurtosis”, The Statistician 47(1), 1998.

KurtosisSample

func (s Series) KurtosisSample(n int) Series

KurtosisSample returns the rolling sample excess kurtosis of s over windows of n bars, G2 = ((n+1) g2 + 6)(n-1)/((n-2)(n-3)), where g2 is the population excess kurtosis of the n values ending at bar t (see Kurtosis). This is the form of Excel’s KURT, pandas’ rolling(n).kurt(), pandas-ta’s kurtosis, scipy’s kurtosis(bias=False) and TA-Lib’s KURTOSIS (2026 main); TradeStation’s Kurtosis is believed to be Excel’s form, but its code is unpublished. On a window that is not flat, the conversion is KurtosisSample = ((n+1) Kurtosis + 6)(n-1)/((n-2)(n-3)), and zero crossings depend on the form: in exact arithmetic, KurtosisSample < 0 exactly when Kurtosis < -6/(n+1). Pearson’s kurtosis m4 / m2^2 (R’s moments::kurtosis, scipy’s fisher=False) is Kurtosis + 3.

KurtosisSample is unitless. Apply it to returns or another stationary series: on price levels it mostly measures the shape of the trend. The sampling error is large on short windows; for normal data the standard error is about sqrt(24/n).

The warm-up reads 0, not NaN: bars 0 to n-2 are 0 and the first valid index is n-1; trim or mask it before using the values as data, as for Kurtosis. The minimum n is 4: an n below 4 or above len(s) gives len(s) zeros, aligned with s. A flat window, whose values are all equal (-0 and +0 count as equal), gives +0, a sentinel rather than a measurement, written directly rather than converted from the 0 that Kurtosis gives; pandas’ rolling kurt gives -3 there (pandas 1.5 and later; older versions give NaN), while scipy and TA-Lib’s KURTOSIS give NaN and Excel #DIV/0!. pandas also gives NaN on varying windows of tiny scale (its threshold is a sum of squared deviations of n*1e-14), where this gives the scale-free value. A NaN or an infinity at bar p gives NaN at bars max(p, n-1) to min(p+n-1, len(s)-1), the windows that hold it, and later windows recover. Finite windows of any magnitude, up to math.MaxFloat64, give finite values.

Each bar is computed from its own window in O(n) time, so the value at a bar does not depend on where a slice of history starts.

Source: Joanes and Gill, “Comparing measures of sample skewness and kurtosis”, The Statistician 47(1), 1998 (the g2 and G2 forms).

MeanAbsDeviation

func (s Series) MeanAbsDeviation(n int) Series

MeanAbsDeviation returns the rolling mean absolute deviation of s over windows of n bars: A/n, where A is the sum of |x - mean| over the n values ending at bar t and mean is the mean of those values (not a moving average or an EMA). It is not the median absolute deviation. This is TA-Lib’s AVGDEV, Tulip’s md, pandas-ta’s mad, TradingView’s ta.dev and Sherry’s “mean deviation” (TA-Lib’s AVGDEV needs n >= 2, and TradingView’s ta.dev skips na values where this gives NaN), and it is the denominator of CCI, which divides by 0.015 times the mean absolute deviation of the typical price.

The result is in the units of the input. In exact arithmetic MeanAbsDeviation(n) <= StdDev(n) <= sqrt(n/2) * MeanAbsDeviation(n). In binary64 the bounds hold to rounding, or to about 1e-16 times sqrt(n) times |mean|/sd relative where the mean dwarfs the spread (StdDev’s running moments lose digits there), while the standard deviation stays below about 1.3e154/sqrt(n) and the deviations from the window mean stay above about 1e-154 in magnitude. Outside that range StdDev’s squared deviations underflow or overflow: Series{1e-170, 3e-170} gives MeanAbsDeviation(2) 1e-170, to rounding, but StdDev(2) 0.

The warm-up reads 0, not NaN: bars 0 to n-2 are 0 and the first valid index is n-1. Trim or mask the warm-up before use: 0 is also the flat-window value. The minimum n is 1: an n below 1 or above len(s) gives len(s) zeros, aligned with s. A flat window, whose values are all equal (-0 and +0 count as equal), gives +0, written directly (here it is also the exact value); at n = 1 every finite window is flat. A NaN or an infinity at bar p gives NaN at bars max(p, n-1) to min(p+n-1, len(s)-1), the windows that hold it, and later windows recover. A finite window gives a finite value unless the exact value is within rounding of math.MaxFloat64, when it can round to +Inf.

Each bar is computed from its own window in O(n) time, so the value at a bar does not depend on where a slice of history starts.

Source: Sherry, “A more conservative estimate of risk”, Stocks & Commodities V7:8 (1989).

RMS

func (s Series) RMS(n int) Series

RMS returns the root-mean-square of n values. Values of n less than two, or greater than the series length, return a zero-filled series of the same length.

The warm-up reads 0, not NaN: bars 0 to n-2 are 0.

Skewness

func (s Series) Skewness(n int) Series

Skewness returns the rolling population skewness of s over windows of n bars: the Fisher-Pearson coefficient g1 = m3 / m2^(3/2) of the n values ending at bar t, where m_k = (1/n) sum (x - mean)^k is their k-th central moment about their own mean. This is the form of Natenberg and Sherry, scipy’s skew(bias=True) and R’s moments::skewness. SkewnessSample is the adjusted form of Excel and pandas, and the conversion between them is SkewnessSample = Skewness * sqrt(n(n-1))/(n-2), so the sign of skewness is the same in every form. Kaufman’s skewness and kurtosis formulas in Stocks & Commodities (September 2026) are non-standard rescalings.

Skewness is unitless. Apply it to returns or another stationary series: on price levels it mostly measures the shape of the trend. The sampling error is large on short windows; for normal data the standard error is about sqrt(6/n). To rounding, |Skewness| <= (n-2)/sqrt(n-1); at n = 2 every window that is not flat gives exactly 0, so the minimum is formal.

The warm-up reads 0, not NaN: bars 0 to n-2 are 0 and the first valid index is n-1. Since 0 is also the skewness of a symmetric window, trim (Window(-(n-1))) or mask (SetN(n-1, math.NaN())) the warm-up before using the values as data. The minimum n is 2: an n below 2 or above len(s) gives len(s) zeros, aligned with s. A flat window, whose values are all equal (-0 and +0 count as equal), gives +0, a sentinel rather than a measurement, as in Standardise and Correlation. A NaN or an infinity at bar p gives NaN at bars max(p, n-1) to min(p+n-1, len(s)-1), the windows that hold it, and later windows recover. Finite windows of any magnitude, up to math.MaxFloat64, give finite values.

Each bar is computed from its own window in O(n) time, so the value at a bar does not depend on where a slice of history starts.

Sources: Natenberg, Option Volatility and Pricing (2nd ed., 2015), Appendix B; Sherry, “A more conservative estimate of risk”, Stocks & Commodities V7:8 (1989); Joanes and Gill, “Comparing measures of sample skewness and kurtosis”, The Statistician 47(1), 1998.

SkewnessSample

func (s Series) SkewnessSample(n int) Series

SkewnessSample returns the rolling sample skewness of s over windows of n bars, the adjusted Fisher-Pearson coefficient G1 = g1 sqrt(n(n-1))/(n-2), where g1 is the population skewness of the n values ending at bar t (see Skewness). This is the form of Excel’s SKEW, pandas’ rolling(n).skew(), pandas-ta’s skew and scipy’s skew(bias=False); TradeStation’s Skew is believed to be Excel’s form, but its code is unpublished. The conversion is SkewnessSample = Skewness * sqrt(n(n-1))/(n-2), so the sign of skewness is the same in both forms.

SkewnessSample is unitless. Apply it to returns or another stationary series: on price levels it mostly measures the shape of the trend. The sampling error is large on short windows; for normal data the standard error is about sqrt(6/n). Carver’s skew forecast, a recipe under “MACD” in the manual, reads it on percentage returns over 60 to 240 days.

The warm-up reads 0, not NaN: bars 0 to n-2 are 0 and the first valid index is n-1; trim or mask it before using the values as data, as for Skewness. The minimum n is 3: an n below 3 or above len(s) gives len(s) zeros, aligned with s. A flat window, whose values are all equal (-0 and +0 count as equal), gives +0, a sentinel rather than a measurement; pandas’ rolling skew also gives 0 there (pandas 1.5 and later; older versions give NaN), while scipy gives NaN and Excel #DIV/0!. pandas also gives NaN on varying windows of tiny scale (its threshold is a sum of squared deviations of n*1e-14), where this gives the scale-free value. A NaN or an infinity at bar p gives NaN at bars max(p, n-1) to min(p+n-1, len(s)-1), the windows that hold it, and later windows recover. Finite windows of any magnitude, up to math.MaxFloat64, give finite values.

Each bar is computed from its own window in O(n) time, so the value at a bar does not depend on where a slice of history starts.

Source: Joanes and Gill, “Comparing measures of sample skewness and kurtosis”, The Statistician 47(1), 1998 (the g1 and G1 forms).

Standardise

func (s Series) Standardise(n int) Series

Standardise returns the z-score of each value against the mean and population standard deviation of its own trailing n-bar window.

For i >= n-1 the result at i is (s[i] - mean)/sigma, where mean and sigma are the mean and the population standard deviation of the window s[i-n+1..i], and the result is 0 when sigma is 0. The computation is ZScore’s rolling-variance sequence with val = s[i] at each bar, so s.Standardise(n)[i] equals s.ZScore(n, s[i])[i] bit for bit, where ZScore scores one constant against every window. Standardise does that work in O(1) amortised time per bar rather than O(len(s)^2) in total.

The first valid index is n-1 and the minimum length is n; earlier outputs are zero. Any n < 2 or n > len(s) gives aligned zeros, a zero-filled series of the same length.

The valid-data class is window-local. A NaN or ±Inf at index p gives NaN for exactly the outputs whose windows contain it, indices max(p, n-1) to min(p+n-1, len(s)-1), and later windows recover. A flat window, in which every value is equal, gives exactly 0. The result is not bitwise start-invariant, because the rolling variance depends on the history before the window. Windows whose sum of squared deviations overflows (deviations above about 1e154/sqrt(n)) give NaN or a signed zero while the extreme values are in the window, and later windows recover.

For the n-1 convention of Ernest Chan’s MATLAB and R listings (his Python listing divides by n) and of pandas-ta’s zscore with talib=False, multiply the result by math.Sqrt(float64(n-1)/float64(n)). The population standard deviation used here is that of TA-Lib’s STDDEV and TradingView’s ta.stdev (biased = true, the default).

Uses include Chan’s pair rules on a spread or ratio (for example long at z <= -2, short at z >= 2 and exit at |z| <= 1) and Markos Katsanos’s Z = (C - MA)/SD. Later relative recipes build on it: the z-score of a log ratio, the %b of a relationship line, rolling co-movement, and the divergence of two z-scores. Jacinta Chan’s BBZ is z := c.Standardise(21) with the state z.GT(1).Sub(z.LT(-1)), and her adjustable-bands ABZ’, whose window varies by bar, is the recipe under StdDevVariable.

Reverre’s ratio mispricing, under “Compositions” in the “Relative: Channels” category, uses it. The later recipes that use it will be listed under “Compositions” in the “Relative: Channels”, “Relative: Ratio and Relative Strength”, “Relative: Regression and Spread” and “Relative: Correlation” categories as they are added.

Verification: a Python reference (exact arithmetic) and an R reference (two-pass population variance) agree, and s.Standardise(n)[i] equals the existing s.ZScore(n, s[i])[i] bit for bit for every i and every n >= 2.

StdDev

func (s Series) StdDev(n int) Series

StdDev returns the population standard deviation of each n values. Values of n less than two, or greater than the series length, return a zero-filled series of the same length.

The divisor is n, as in TA-Lib’s STDDEV, TradingView’s ta.stdev and Excel’s STDEV.P. pandas, pandas-ta’s stdev with talib=False and Excel’s STDEV.S divide by n-1: that is StdDevSample(n), equal to StdDev(n) times sqrt(n/(n-1)) to rounding on ordinary prices; see StdDevSample. UpsideSemiDeviation and DownsideSemiDeviation split the deviation by side of the window’s mean: N+ sigma+^2 + N- sigma-^2 = n StdDev(n)^2 to rounding, with N+ and N- the counts above and below it.

The warm-up reads 0, not NaN: bars 0 to n-2 are 0, which reads as no volatility.

Precision: StdDev updates running moments as values enter and leave the window, and their rounding error scales with the level of the values, not with their spread. The relative error is up to about 1e-16 * sqrt(n) * |mean| / sd, with mean and sd those of the window, or a few units in the last place where the mean is small. On ordinary prices that is far below anything a trading rule would see, but it grows with the level: at n = 10, on prices that move by about 0.01 a bar, the worst error measured is about 1e-12 near 100, 2e-10 near 1e4, 1e-8 near 1e6 and 1e-6 near 1e8. StdDevSample, StdDevVariable and MeanAbsDeviation centre and scale each window first and do not have this error. The methods built on StdDev or its running moments inherit it: Bollinger, ZScore, Standardise, VIDYA and RelativeVolatilityIndex among them. The kernel is kept as it is, since a fix would change those methods’ bits (GLE-472, 10 October 2026).

StdDevSample

func (s Series) StdDevSample(n int) Series

StdDevSample is the sample standard deviation of the last n values, with divisor n - 1: Excel’s STDEV.S, pandas’s rolling(n).std(), pandas-ta’s stdev with talib=False and Julia’s std. StdDev is the population form, with divisor n. StdDevSample(n) and StdDev(n) times sqrt(n / (n - 1)) agree to rounding on ordinary data, and to about 1e-16 times sqrt(n) times |mean|/sd relative in general: StdDev’s running moments lose digits when the mean dwarfs the spread, overflow once the standard deviation passes about 1.3e154/sqrt(n), and lose precision once the deviations from the window mean fall below about 1e-154 in magnitude, reaching 0 near 1e-162, where this kernel stays accurate. The two use different kernels, so no bit-for-bit identity holds between them. XStdDev is the exponentially weighted form, and StdDevVariable the population form with its period read from a series at each bar.

The warm-up reads 0, not NaN: bars 0 to n-2 are 0 and the first valid index is n-1. Trim or mask the warm-up before use, since 0 is also a legitimate reading. n < 2 or n > len(s) gives len(s) zeros.

NaN footprint: a window holding a NaN or an infinity gives NaN, so a bad value at p makes bars max(p, n-1) to min(p+n-1, len(s)-1) NaN; later windows recover. This test comes before the flat test, so a window of equal infinities is NaN. A finite window whose values are all equal gives +0 exactly. The window is centred and scaled before squaring, so a finite window gives a finite value unless the sample standard deviation itself exceeds math.MaxFloat64, which gives +Inf. Each value depends only on its own window, so the result is bitwise window-local.

StdDevSample(n).Squared() is the sample variance. The close-to-close volatility is the sample standard deviation of the log returns ln(c[t] / c[t-1]), r := c.LogRatio(c.Shift(1)), over the window; multiply by the square root of the number of bars in a year to annualise it. The log return at bar 0 is undefined: Shift pads bar 0 with 0, so r is NaN there and the first window is NaN; take the returns from bar 1, by trimming bar 0 or by setting it with SetN(1, 0) and reading the result from bar n, the first window without it. For TTR’s close estimator over n prices, see YangZhangVolatility.

The cost is O(n) per bar; only the result and one scratch slice of length n are allocated.

UpsideSemiDeviation

func (s Series) UpsideSemiDeviation(n int) Series

UpsideSemiDeviation returns the upside semi-deviation of s over a rolling window of n bars, the root mean square deviation from the window’s mean of the values above that mean. For the n values ending at bar t, with mu their mean, it is sigma+ = sqrt(S+ / N+) (Poster’s equation 1), where S+ is the sum of (x - mu)^2 over the values x > mu and N+ is their number. The inequality is strict: a value equal to mu is on neither side, so N+ + N- <= n, where N- counts the values below mu. The divisor is the side’s own count N+, not n, so sigma+ is the root mean square deviation of the values above the mean. Poster’s indicator deviation bands replace Bollinger’s single standard deviation with one deviation per side of the mean, sigma+ above it and sigma- below it (see DownsideSemiDeviation).

The input is any series, typically an indicator: Poster applies the bands to a fast SMA of the close, to the RSI and to the stochastic.

The mean is the method’s own, computed as in StdDevSample: the window is anchored at its first value, scaled by a power of two and averaged with a compensated sum. It agrees with Average(n) to rounding, not bit for bit, and the method does not call Average. A value counts as above the mean when its scaled centred deviation is positive and as below it when the deviation is negative, by strict comparisons with 0 rather than by math.Signbit, so a deviation of -0, like one of +0, is on neither side. Membership agrees with exact arithmetic on windows of integers, or of multiples of a common power of two, whose span, counted in units of that power and multiplied by n, is below 2^52: there the differences from the first value and their sum are exact, and the rounding of the mean does not move any value to another side. Ties with a decimal mean are fragile: a value equal to the decimal mean may count on either side or on neither, on data with few decimals such as an RSI to two decimals or prices to the tick, because binary floating point cannot hold the decimal mean and no mean of the stored values, whether this one, Average(n), a naive sum or their exact mean, finds such ties every time. That can move N+ by one and sigma+ by a step, a property of the count normalisation, not of this method.

In exact arithmetic, other conventions differ from this one by a factor that changes from bar to bar with N+, not by a constant: the Markowitz semi-deviation divides by n and is sigma+ times sqrt(N+/n), and the convention of Kahn and Stefek divides by n-1, which gives sigma+ times sqrt(N+/(n-1)). The two sides add up to the variance, to rounding: N+ sigma+^2 + N- sigma-^2 = n StdDev(n)^2 = (n-1) StdDevSample(n)^2, with a value at the mean adding to neither side. On a window symmetric about its mean with no value at the mean the two sides are equal, bit for bit in the region above where membership agrees with exact arithmetic, and both equal StdDev(n) to rounding, so bands built on the two semi-deviations about Average(n) are then Bollinger’s to rounding. Negation swaps the sides exactly: s.Mul(-1).DownsideSemiDeviation(n) equals s.UpsideSemiDeviation(n) bit for bit, for every input. At n = 2 both sides are |s[t] - s[t-1]|/2 for t >= 1 while the difference does not overflow, and StdDev(2) bit for bit while that difference lies between about 3e-154 and 1.9e154 in magnitude, outside which StdDev’s running moments underflow or overflow; a window whose difference overflows, such as Series{math.MaxFloat64, -math.MaxFloat64}, is scaled first here and gives math.MaxFloat64 where StdDev(2) gives +Inf.

Poster’s bands are four lines about the simple moving average of the indicator, two above it at k1 and k2 times sigma+ and two below it at k1 and k2 times sigma-, with k2 > k1. For an indicator x and a band period n,

mid := x.Average(n)
sp, sn := x.UpsideSemiDeviation(n), x.DownsideSemiDeviation(n)
obLwr, obUpr := mid.Add(sp.Mul(k1)), mid.Add(sp.Mul(k2))
osUpr, osLwr := mid.Sub(sn.Mul(k1)), mid.Sub(sn.Mul(k2))

gives the overbought zone from obLwr to obUpr and the oversold zone from osLwr to osUpr. The midline and both semi-deviations take the same n, as Poster’s mu is the average of the indicator over the band period; the indicator’s own period and any smoothing of it, which Poster optimises separately, are the caller’s. Another midline, an EMA for instance, can replace mid, but sigma+ and sigma- still measure deviations from the window’s simple mean. The factors are typically 1.5 to 2.5, and the dual Bollinger Bands of Poster’s Figure 1 use 1.5 and 2.0. On bars 0 to n-2, Average passes the input through and both semi-deviations are 0, so all four bands equal the input there, as Bollinger’s do. Nothing clips the bands, so on a bounded oscillator such as the RSI or the stochastic they can leave its range of 0 to 100. Poster’s trading rules (entry into a zone, a limit on the indicator’s change, stops and targets, and a filter on the width of Bollinger Bands on price) are not part of the recipe.

The result lies between 0 and the largest |x - mu| of the window, to rounding. A finite window gives a finite value unless the exact value reaches math.MaxFloat64 plus half a unit in its last place, where binary64 rounds to +Inf, or comes within rounding of that, when it is +Inf: with n = len(s), Series{math.MaxFloat64, -math.MaxFloat64, 1}, whose exact downside value is math.MaxFloat64 plus 1/3, gives math.MaxFloat64. With n = len(s), Series{math.MaxFloat64, -math.MaxFloat64} gives math.MaxFloat64 on both sides, as a window whose differences overflow is scaled first. Series{math.MaxFloat64, -math.MaxFloat64, math.MaxFloat64}, with n = len(s), gives 1.1984620899082105e308 here, 2/3 of math.MaxFloat64, and +Inf from DownsideSemiDeviation, where the exact value is 4/3 of math.MaxFloat64. At the other end, a window that is not flat gives 0 only when the exact value is about 2.5e-324 or less, and with n = len(s), Series{0, 5e-324} gives 0 on both sides.

XBeta

func (s Series) XBeta(s2 Series, n int) Series

XBeta returns the exponentially weighted (EW) beta of s on s2, for the span n, so that alpha = 2 / (n + 1): s2 is the regressor and s the response. The beta is the EW covariance of the two series divided by the EW variance of s2.

Series are paired by position only; see “Time alignment” in the manual. The series are right-aligned at the most recent bar and the beta is taken over the overlap. The result is as long as s, and bars before the overlap are 0.

The first n bars of the overlap are the seed: the state starts as the population means, the variance of s2 and the covariance of that window, exactly as XCorrelation seeds its state, and the seed bar is the least-squares slope of that window, as Beta gives there. Each later bar updates the state as in West (1979). The weights are alpha(1 - alpha)^k on the bars after the seed, k being the bar’s distance back from t, and the seed window shares the weight (1 - alpha)^(t-n+1) equally, so the seed’s weight decays by the factor 1 - alpha per bar. The value at bar t includes bar t.

Flat rule: the value is 0 wherever the trailing n values of either input in the overlap, up to and including the current bar, are all equal (+0 and -0 count as equal), which is where XCorrelation is 0 by flatness; a NaN bar gives NaN. The rule changes the output only, and the state goes on updating. The value is also 0 where the EW variance of s2 is not positive. A nearly flat regressor makes beta explode, as with Beta, and a quotient that overflows is left as +Inf or -Inf.

The warm-up reads 0, not NaN: bars 0 to n-2 are 0 and the first valid index is n-1, counted from the start of the overlap; bars before the overlap are 0. Trim or mask the warm-up before use, since 0 is also a legitimate reading. Invalid parameters give a series of zeros the same length as s: n < 2, an empty overlap, or n greater than the overlap’s length.

A NaN or infinity in the seed window of either series, or a seed variance of s2 or covariance that overflows, makes every value from the seed bar on NaN. Later, a NaN or infinite input in either series, or a mean, variance of s2 or covariance that overflows, makes the value at that bar and every later bar NaN, without recovery. The method is recursive, so under a sliding slice of history the value depends on where the slice starts.

Beta is unitless only if s and s2 share units; otherwise it is in the units of s per unit of s2. See XVariance for the conventions and sources.

XCovariance

func (s Series) XCovariance(s2 Series, n int) Series

XCovariance returns the exponentially weighted (EW) population covariance of s and s2 about their EW means, for the span n, so that alpha = 2 / (n + 1). s.XCovariance(s, n) equals s.XVariance(n), and over the overlap the covariance is symmetric in the two series.

Series are paired by position only; see “Time alignment” in the manual. The series are right-aligned at the most recent bar and the covariance is taken over the overlap. The result is as long as s, and bars before the overlap are 0.

The first n bars of the overlap are the seed: the state starts as the population means and covariance of that window, exactly as XCorrelation seeds its state. Each later bar updates the state as in West (1979), with the product of the two deviations in the covariance. The weights are alpha(1 - alpha)^k on the bars after the seed, k being the bar’s distance back from t, and the seed window shares the weight (1 - alpha)^(t-n+1) equally, so the seed’s weight decays by the factor 1 - alpha per bar. The value at bar t includes bar t.

A NaN or infinity in the seed window of either series, or a seed covariance that overflows, makes every value from the seed bar on NaN. Later, a NaN or infinite input in either series, or a mean or covariance that overflows, makes the value at that bar and every later bar NaN, without recovery. The variance of either series plays no part, so the result is symmetric in s and s2 even where one of them has a variance beyond the float64 range. The method is recursive, so under a sliding slice of history the value depends on where the slice starts.

The covariance is in the product of the two inputs’ units. See XVariance for the conventions and sources.

XStdDev

func (s Series) XStdDev(n int) Series

XStdDev returns the exponentially weighted (EW) standard deviation of s: the square root of XVariance(n), for the span n, so that alpha = 2 / (n + 1). It has the same seed, the population moments of the first n-bar window as in XCorrelation, and the same update (West, 1979), so the weights, the non-finite rule and the recursion are those of XVariance, and the seed’s weight decays by the factor 1 - alpha per bar. The value at bar t includes bar t.

The warm-up reads 0, not NaN: bars 0 to n-2 are 0 and the first valid index is n-1. Trim or mask the warm-up before use: a zero standard deviation reads as no volatility. Invalid parameters give a series of zeros the same length as s: n < 2, n > len(s), or an empty s.

A NaN or infinity in the seed window, or a seed window whose variance overflows, makes every value from bar n-1 on NaN. Later, a NaN or infinite input, or a mean or variance that overflows, makes the value at that bar and every later bar NaN, without recovery. The method is recursive, so under a sliding slice of history the value depends on where the slice starts.

The result is in the units of s. Carver’s span is 32 (35 is also used) and he annualises the standard deviation of daily returns by 16. RiskMetrics’ decay factor 0.94 corresponds to alpha = 0.06 in XStdDevAlpha. As a forecast for the next bar (RiskMetrics’ convention) the value at bar t needs no shift, while forecasters who index it at t + 1 shift it by one bar. See XVariance for the sources and the pandas conventions. For the equally weighted sample standard deviation of a trailing n-bar window, see StdDevSample.

XStdDevAlpha

func (s Series) XStdDevAlpha(alpha float64) Series

XStdDevAlpha returns the exponentially weighted (EW) standard deviation of s: the square root of XVarianceAlpha(alpha), for the smoothing factor alpha. It is seeded with the first value and updated as in West (1979), so the weights, the non-finite rule and the recursion are those of XVarianceAlpha, and the seed’s weight decays by the factor 1 - alpha per bar. The value at bar t includes bar t.

The alpha form has no warm-up: every bar is valid and bar 0 is 0. alpha must be finite with 0 < alpha <= 1; otherwise the result is zeros the same length as s. An empty s gives an empty result.

A NaN or infinite s[0] makes every value NaN. Later, a NaN or infinite input, or a mean or variance that overflows, makes the value at that bar and every later bar NaN, without recovery. The method is recursive, so under a sliding slice of history the value depends on where the slice starts.

The result is in the units of s. RiskMetrics’ decay factor 0.94 is alpha = 0.06, and Carver annualises the standard deviation of daily returns by 16. See XVariance and XVarianceAlpha for the conventions and sources.

XVariance

func (s Series) XVariance(n int) Series

XVariance returns the exponentially weighted (EW) population variance of s about its EW mean, for the span n, so that alpha = 2 / (n + 1). It is the EW variance of Carver, and RiskMetrics’ weighting applied to deviations from an EW mean (RiskMetrics itself uses a zero mean, below), with the mean and the variance updated as in West (1979).

The first n bars are the seed: the state starts as the population mean and variance of that window, exactly as XCorrelation seeds its state. Each later bar t takes d = s[t] - mean, adds alphad to the mean and sets the variance to (1 - alpha)(variance + alphadd). The weights are alpha(1 - alpha)^k on the bars after the seed, k being the bar’s distance back from t, and the seed window shares the weight (1 - alpha)^(t-n+1) equally, so the seed’s weight decays by the factor 1 - alpha per bar. The value at bar t includes bar t.

The warm-up reads 0, not NaN: bars 0 to n-2 are 0 and the first valid index is n-1. Trim or mask the warm-up before use: a zero variance reads as no volatility. Invalid parameters give a series of zeros the same length as s: n < 2, n > len(s), or an empty s.

The variance is in the units of s squared.

Conventions: Carver uses a span of 32 (also 35) on daily returns and annualises by multiplying the standard deviation by 16 (so the variance by 256). RiskMetrics’ decay factor 0.94 corresponds to alpha = 0.06 in XVarianceAlpha, and the zero-mean RiskMetrics variance is r.Mul(r).XAverageAlpha(0.06). The sample (bias-corrected) variance is n/(n-1) times this value in exact arithmetic. As a forecast for the next bar (RiskMetrics’ convention) the value at bar t needs no shift, while forecasters who index it at t + 1 shift it by one bar. pandas ewm(adjust=False) is the alpha form; pandas’ default adjust=True, used by López de Prado’s getDailyVol and pysystemtrade, weights the first bars differently and is not reproduced, and their default sample correction (bias=False) keeps their standard deviation about sqrt(n/(n-1)) above XStdDev.

Sources: RiskMetrics Technical Document (J.P. Morgan/Reuters, 1996), section 5.2; West, “Updating mean and variance estimates: an improved method”, Communications of the ACM 22(9), 1979; Carver, Advanced Futures Trading Strategies (2023); and, for decayed sums and online moments, Isichenko, Quantitative Portfolio Management (2021).

XVarianceAlpha

func (s Series) XVarianceAlpha(alpha float64) Series

XVarianceAlpha returns the exponentially weighted (EW) population variance of s about its EW mean, for the smoothing factor alpha, seeded with the first value: the mean is s[0] and the variance is 0 at bar 0. It equals pandas Series.ewm(alpha=alpha, adjust=False).var(bias=True) to rounding. XVariance seeds with the first n-bar window instead, so XVariance(n) is not XVarianceAlpha(2 / (n + 1)).

Each later bar t takes d = s[t] - mean, adds alphad to the mean and sets the variance to (1 - alpha)(variance + alphadd), as in West (1979). The weight of bar t-k is alpha(1 - alpha)^k for k < t, and bar 0 carries the remaining (1 - alpha)^t, so the seed’s weight decays by the factor 1 - alpha per bar. The value at bar t includes bar t.

The alpha form has no warm-up: every bar is valid and bar 0 is 0. alpha must be finite with 0 < alpha <= 1; otherwise the result is zeros the same length as s. An empty s gives an empty result.

The variance is in the units of s squared. Under this first-value seed the sample (bias-corrected) variance is 1/(1 - sum of the squared weights) times this value, a factor that varies by bar, unlike XVariance’s n/(n-1).

RiskMetrics’ decay factor lambda = 0.94 is alpha = 0.06: pass 0.06. A constant expression 1 - 0.94 is exact in Go and also gives 0.06, but 1 - lambda formed at run time from a float64 variable lambda = 0.94 is 0.060000000000000053. The zero-mean RiskMetrics variance is r.Mul(r).XAverageAlpha(0.06); GARCHVariance is its GARCH(1,1) generalisation with a long-run variance term. As a forecast for the next bar (RiskMetrics’ convention) the value at bar t needs no shift, while forecasters who index it at t + 1 shift it by one bar. pandas’ default adjust=True, used by López de Prado’s getDailyVol and pysystemtrade, weights the first bars differently and is not reproduced. See XVariance for the sources.

ZScore

func (s Series) ZScore(avgLen int, val float64) Series

ZScore returns the z-score of the constant val against the mean and population standard deviation of each trailing avgLen-bar window.

The z-score is a statistical measure of how many standard deviations a value is from the mean. Windows with zero standard deviation, and avgLen values outside the bounds of the series, yield zero. For the z-score of each value against its own window, use Standardise.

The output is in standard deviations. The warm-up reads 0, not NaN: bars 0 to avgLen-2 are 0, which reads as at the mean.

Percent Rank

PercentRank

func (s Series) PercentRank(val float64, n int) Series

PercentRank computes the percentile rank of a value within a rolling window of n bars. Returns a value between 0 and 1 indicating where the given value falls within the sorted window, or -1 if the value is outside the window’s range. Consistent with TradeStation’s PercentRank function.

The output is on a 0 to 1 scale, despite the name, and is -1 when val lies outside the window’s range. The warm-up reads 0, not NaN: bars 0 to n-2 are exactly 0, the first valid index is n-1 and the minimum length is n. 0 reads as the bottom of the window; trim (Window(-(n-1))) or mask (SetN(n-1, math.NaN())) the warm-up first.

R-Squared

R2

func (s Series) R2(n int) Series

R2 computes a least-squares regression on the last <n>{=html} values at each point in the series. Where there are insufficient values, the result is zero. The underlying calculation uses numerically stable, centred coordinates so R2 remains accurate for windows with a large common offset (for example, prices near 1e9). Finite constant and singleton windows yield zero. A window containing NaN or infinity yields NaN; finite windows recover.

The output is on a 0 to 1 scale, not a percentage. Its warm-up reads 0, not NaN, on bars 0 to n-2, and 0 reads as no linear fit.

Spearman Rank Correlation

SpearmanRankCorrelation

func (s Series) SpearmanRankCorrelation(n int, s2 ...Series) Series

SpearmanRankCorrelation calculates the Spearman Rank Coefficient (correlation) over a rolling window of n bars. If a second series is passed in as a parameter, the receiver is ranked against it; only the first such series is used and any further ones are ignored. If no second series is provided, each window of the receiver is ranked against a sorted (low to high) copy of itself, which provides an indication of the strength of a trend (since the series is being ranked against itself in ascending order). A second series is right-aligned with the receiver at the most recent bar: with m = min(len(s), len(s2[0])), the last m bars of each are paired, the correlation of each pair of windows is written at the position of the window’s latest receiver bar, and values of either series before the overlap are never read. The result always has the receiver’s length. It is zero before the first complete window, and entirely zero when n < 2 or n exceeds the bars available (m with a second series, including an empty one, or len(s) without). A window in which any value of either series is NaN, +Inf or -Inf gives NaN at that position. Tied values receive the average of their ranks, and a window whose ranks are all tied gives 0. The rank scratch is allocated once per call, so the number of allocations does not grow with len(s) or n. Series are paired by position only; see “Time alignment” in the manual.

Linear Regression

LinRegValue

func (s Series) LinRegValue(n int) Series

LinRegValue returns the value at the current bar of the least-squares line fitted to the trailing n values.

For i >= n-1 the line is the ordinary least-squares fit of the window s[i-n+1:i+1] against bar number, evaluated at the window’s last bar, the current one. The line passes through the window’s mean at the middle of the window, (n-1)/2 bars before the current bar, so the output is the n-bar Average plus the n-bar Slope (in units of s per bar) times (n-1)/2: computed as s.Average(n)[i] + s.Slope(n)[i]h, with h = float64(n-1)/2 and the product rounded to binary64 before the addition. Hence the relation LinRegValue(n) - Average(n) = Slope(n)(n-1)/2 holds, to rounding, from the first valid index on. A flat window gives Average’s value of that window, its Slope being exactly 0, and an exactly linear window gives the current bar’s value, to rounding, and exactly on a small dyadic ramp such as 128 + 0.25*i.

The first valid index is n - 1 and the minimum length is n; every output before the first valid index is 0 (positive zero). An empty input, n < 2 or n > len(s) (math.MaxInt and math.MinInt included) returns a fresh all-zero series of len(s) in O(len(s)) time, allocating only the result and without computing any average or slope.

A NaN, +Inf or -Inf at index p gives NaN on [max(p, n-1), p+n-1], because Slope is NaN for any window that holds one; later outputs recover. Finite inputs whose window sum overflows binary64 can give a non-finite output, whatever the composition of Average and Slope gives; this is not special-cased.

It is not start-invariant: Average keeps a rolling compensated sum, so dropping bars from the front of s can change later outputs in the last bits, although in exact arithmetic each output depends only on the n values of its own window. Valid data: any finite series; large common offsets are handled, because Slope fits centred coordinates and Average sums with compensation. It costs O(n) per bar, through Slope, and O(len(s)*n) per call. It only reads s, so it is safe for concurrent use on a shared input, and its result is freshly allocated, sharing no storage with s. For valid n it reuses Slope’s result and allocates nothing beyond one Average call and one Slope call.

It computes MetaStock’s LinearReg(x, n), the end-point value of the n-bar linear regression line. Deviation from MetaStock: every output before the first valid index is 0, where MetaStock has no value. The least-squares moving average (LSMA) is the same quantity under another name, and TA-Lib’s LINEARREG computes it too (not bit for bit: TA-Lib sums without compensation and may fuse the multiply and add). The time-series forecast, one bar ahead on the same line, is s.LinRegValue(n).Add(s.Slope(n)); TA-Lib’s TSF is this quantity. PolyFit(n, 1, 0) is this end point and PolyFit(n, 1, 1) that forecast, the line as the least-squares polynomial of degree 1, both to rounding (see ExampleSeries_PolyFit_tsf). In exact arithmetic LinRegValue(n) is also 3WAverage(n) - 2Average(n), the linearly weighted average’s known link to the regression end point, which holds to rounding in binary64. Ehlers’s projected moving average, PMA(n), adds half a bar more of slope (Average + Slopen/2), so from the first valid index on LinRegValue(n) equals PMA(n) - 0.5Slope(n) in exact arithmetic, and to rounding in binary64.

Markos Katsanos, Intermarket Trading Strategies (Wiley, 2008), uses it in his relative-strength trend system. His text (page 184) reads “Buy when the RS line crosses over its linear regression line of six days ago … the 200-day regression line was the most profitable for long entry and the 50-day for long exit”, and his code (Appendix A, pages 304-305) reads “RS:=MOV((C/SEC2),3,E); LR1:=REF(LinearReg(RS,D1+D2,E,3),-D2);”, with D1 = 200 and D2 = 15 for the buy and D1 = 50 and D2 = 6 for the sell. On a trimmed relative-strength line rs (his RS), his LR1 translates to rs.LinRegValue(215).Shift(15) for the buy and rs.LinRegValue(56).Shift(6) for the sell. The book’s four-argument “LinearReg(RS,D1+D2,E,3)” does not match MetaStock’s two-argument function and is unexplained; the translation reads it as the two-argument “LinearReg(RS, D1+D2)”. The text’s regression line “of six days ago” matches the sell code’s lag (D2 = 6), while the buy code lags 15 bars, and the text’s 200-day and 50-day lines are D1, whereas the code fits D1+D2 bars (215 and 56).

The warm-up 0 before the first valid index reads as a regression line at zero. On a positive series, Katsanos’s level rule “RS > (1+3/100)*LR1” therefore holds on every warm-up bar, and his exit cross “CROSS(LR1,RS)”, the line crossing above the series, can fire falsely at the first valid bar, where the line jumps from 0 to a value above the series. (A cross of a positive series above the line cannot fire on the warm-up, because the series is already above the zero line.) Trim (Window(-(n-1))) or mask (SetN(n-1, math.NaN())) the warm-up before such rules. The Shift(15) of his buy form delays the first valid value to index n-1+15 (the Shift(6) of the sell form to n-1+6), and Shift fills its start with 0, so trim or mask n-1+15 (or n-1+6) bars of the shifted line.

Verification: two independent references agree with it to rounding, a Python reference in exact rational arithmetic and a Julia reference that solves each window’s least-squares problem with Julia’s own solver (the LinearAlgebra standard library) in BigFloat. From the first valid index on it equals s.Average(n).Add(s.Slope(n).Mul(float64(n-1) / 2)), the composition of Average and Slope, bit for bit, NaN matching NaN; and the negation and power-of-two scaling identities hold bit for bit, apart from the sign of a zero and of a NaN: s.Mul(-1).LinRegValue(n) equals s.LinRegValue(n).Mul(-1), and a power-of-two scale of s that neither overflows nor underflows scales the result by the same power of two.

PolyFit

func (s Series) PolyFit(n, degree int, ahead float64) Series

PolyFit returns the value of the least-squares polynomial of the given degree fitted to the last n values, evaluated ahead bars past the newest of them.

At each bar t from n-1 on, it fits the polynomial of degree degree by ordinary least squares to the window s[t-n+1] to s[t], placed at the abscissae 0 to n-1 with the current bar at n-1, and evaluates it at n-1+ahead. The fitted value is linear in the window, and the result is s.FIRSum(PolyFitWeights(n, degree, ahead)), bit for bit: the weights are the exact least-squares weights, each correctly rounded to float64, and each bar is their dot product with its window, each product rounded on its own and summed oldest first from +0, plus 0, so a zero result is +0 and every platform gives the same bits. The origin of the axis does not matter: Meyers’s fit to t = 1 to T, evaluated at T+ahead, is the same number. With x_j the window’s values, w_j the exact weights, S the sum of |w_j x_j| and u = 2^-53, each valid bar lies within (n+2)u S of the exact least-squares value when no weight and no product is subnormal. On bars where no weight and no product is subnormal, the error measured on the reference fixtures is at most 12.3u S, and 13.1u S on a 3,000-bar log walk at n = 200; these are properties of those series, not bounds, and the error grows with n, so only (n+2)u S is guaranteed. No level is subtracted from the window first, unlike Slope(n), so a large common offset costs accuracy relative to the fluctuations about it.

Any finite float64 ahead is accepted, large ones included, though from degree 2 a huge ahead overflows the weights and gives zeros (see PolyFitWeights). Ahead 0 gives the end point of the fitted curve and ahead 1 the one-bar forecast, and at degree 1 ahead 0.5 gives the point PMA(n) projects to. A negative ahead gives a fitted value inside the window, a lagged smoother stamped at the current bar, which is causal and does not repaint: ahead -(n-1) gives the oldest point and ahead -(n-1)/2 the centre, which for odd n is the centred Savitzky-Golay smoother delayed by (n-1)/2 bars. The degree lies between 0 and 4, with n at least degree+1. With n = degree+1 the polynomial interpolates the window, so ahead 0 gives the current value, and degree 0 is the moving average, Average(n) to rounding, for every ahead. The recursive counterpart, a polynomial fitted to the whole past with weights that decay geometrically rather than equally over a window, is the fading-memory family, FadingMemoryLevel and its siblings.

The method is Meyers’s, from “The British Pound, Cubed” (Stocks & Commodities V16:11, November 1998), whose forecast curve P_f, the cubic fitted to the last T closes and extended one bar, is c.PolyFit(T, 3, 1), with T = 40 to 44 in his optimisations. In “IBM, Cubed” (Stocks & Commodities V17:8, August 1999) he fits the midpoint (H+L)/2, as in h.Add(l).Div(2).PolyFit(34, 3, 1), and lets the optimiser pick the degree, 2 or 3, and in “The British Pound Cubed, Redux” (Stocks & Commodities V31:1, January 2013) he takes T from 20 to 60. He solves the raw-power normal equations, which are badly conditioned in floating point, where these weights are exact before their one rounding.

The 1998 article’s Excel sidebar, “Using Excel for a cubic polynomial” by Thom Hartle, fits a cubic to ten days: its fitted column is the end point PolyFit(10, 3, 0), not the forecast, and prints 6.03 on the tenth day of the series 1, 2, 3, 4, 5, 4, 3, 4, 5, 6, where the forecast PolyFit(10, 3, 1) is 41/5. Its chart’s cubic is y = 0.0342x^3 - 0.5944x^2 + 3.3411x - 2.1333 and its line y = 0.4061x + 1.4667, whose slope, 67/165, is Slope(n) at n = 10. Meyers’s percent-reversal and two-day jump rules are trading recipes on the curve, not part of the method; the percent-reversal state is DirectionalChange, applied to the curve, and the rules themselves are not built here. His tests are walk-forward, by the author alone, each on one market and before costs.

To rounding, on finite data: degree 0 matches Average(n) from bar n-1; degree 1 matches Average(n) + Slope(n)*((n-1)/2 + ahead) from bar n-1; degree 1 at ahead 0 matches LinRegValue(n) over the whole series; at ahead 1 it matches the time-series forecast, TA-Lib’s TSF, s.LinRegValue(n).Add(s.Slope(n)), the curve of Meyers’s earlier linear system, over the whole series for n >= 2; at ahead 0.5 it matches PMA(n) from bar n-1; and a polynomial q of degree at most degree, sampled as s[i] = q(i), gives q(t+ahead) at bar t. Bit for bit: n = 1 (degree 0) gives each value plus 0; n = degree+1 at ahead 0, whose weights are 0, …, 0, 1, gives each value plus 0 from bar n-1 on finite data; from bar n-1, degree 0 with n a power of two equals Average(n) on small integers; negating s negates the result, a zero result being +0 either way; and scaling s by a power of two scales the result exactly, as long as nothing overflows, underflows or becomes subnormal. The derivative of the fitted polynomial is not provided; at degree 1 the slope is Slope(n). A log-linear or other transformed fit takes a transformed input, as in c.Log().PolyFit(n, degree, ahead).

The warm-up reads 0, not NaN: bars 0 to n-2 are 0 and the first valid index is n-1. Trim it with Window(-(n-1)) for n >= 2, since Window(0) is empty, or mask it with SetN(n-1, math.NaN()). The minimum length is n, and n = len(s) gives one valid bar, the last. The result is len(s) zeros when n < 1, n > len(s), the degree is outside 0 to 4, n < degree+1, ahead is NaN or infinite or a weight overflows; n > len(s) is checked before any weights are built.

Non-finite values follow IEEE arithmetic and stay in the windows that hold them: a NaN gives NaN, and an infinity gives an infinity with the sign of its weighted term, or NaN when the window’s weighted terms are infinities of both signs or the infinity meets a weight of +0 or -0, since 0*Inf is NaN. Some weights are exactly 0, as at n = degree+1 and ahead 0, so Series{+Inf, 1, 2}.PolyFit(2, 1, 0) is {0, NaN, 2}. Later windows recover. Finite values can overflow: Series{-1e308, 1e308}.PolyFit(2, 1, 1) is {0, +Inf}. There is no flat-window rule: a flat window of x gives x to within the rounding of the weighted terms, (n+2)u S with S = |x| times the sum of the |w_j|. The error is tens of ulps of x at most for n up to 60 and an ahead from -(n-1) to 1 (31 measured), but grows with n (137 ulps of 1e-10 at n = 1000 and degree 0) and with extrapolation (640 ulps of -7.5 at n = 52, degree 4 and ahead 52), and far more for a huge ahead, whose large weights of both signs cancel (a flat window of 1 at degree 2 and ahead 1e10 gives 8192, within (n+2)u S).

The white-noise gain, the square root of the sum of the squared weights, is the factor by which the fit scales independent noise. At degree 3 and ahead 1 it is 1.09 at n = 20, 0.77 at n = 34 and 0.66 at n = 44, while the end point at n = 44 has 0.56 and the 44-bar average 0.15, so the one-bar forecast of a short window amplifies noise rather than smoothing it.

Each valid bar costs O(n) and depends only on its own window, so the result is window-local and start-invariant, bit for bit: dropping k bars from the front leaves every bar from k+n-1 on unchanged, and s[:j] gives the first j values of the result. The weights are built once per call (see PolyFitWeights), so a call costs O(len(s)*n) after them. PolyFit is safe for concurrent use, does not modify s and returns a new series.

Verification: two clean-room references, Python with exact fractions and Julia with exact rationals, compute the exact weights and window fits, confirm that the weights are the correctly rounded exact weights and that every output lies within a measured bound of the exact fit, and their binary64 transcriptions agree with PolyFit bit for bit. NumPy’s polyfit and R’s lm with a raw polynomial are tolerance cross-checks, and from bar n-1 the degree-0 fit of a power-of-two window equals Average(n) bit for bit on small integers.

PolyFitWeights

func PolyFitWeights(n, degree int, ahead float64) []float64

PolyFitWeights returns the n weights, oldest first, of the least-squares polynomial of the given degree fitted to a window of n values and evaluated ahead bars past its newest value.

With the window at the abscissae j = 0 to n-1, X the n by degree+1 matrix whose row j is (1, j, …, j^degree), h = n-1+ahead and e(h) = (1, h, …, h^degree), the weights are, in normal-equation form and in Gram form,

w = e(h)^T (X^T X)^-1 X^T,
w[j] = sum over k = 0 to degree of P_k(u_j) P_k(u_h) / ||P_k||^2,

where u_j = j - (n-1)/2 and u_h = h - (n-1)/2 are the centred abscissae, the squared norm ||P_k||^2 is the sum over j of P_k(u_j)^2, and P_0(u) = 1, P_1(u) = u and P_{k+1}(u) = u P_k(u) - (k^2 (n^2 - k^2) / (4 (4k^2 - 1))) P_{k-1}(u) are the monic discrete orthogonal (Gram) polynomials of the n points. The fitted value at h is w[0]*x[0] + … + w[n-1]*x[n-1] for the window x[0] to x[n-1], so s.FIRSum(PolyFitWeights(n, degree, ahead)) is s.PolyFit(n, degree, ahead).

The weights are computed exactly in rational arithmetic with math/big, from the exact value of ahead, and each is rounded once to the nearest float64, ties to even, by big.Rat’s Rat.Float64, so they are fixed by mathematics rather than by an order of operations and are the same on every platform. A weight that is exactly 0 is +0, and a negative weight that rounds to zero is -0: PolyFitWeights(5, 1, 0) is [-0.2 0 0.2 0.4 0.6], whose second weight is exactly 0, while the second weight of PolyFitWeights(5, 1, 5e-324) is negative and rounds to -0. In exact arithmetic the weights sum to 1 and reproduce every polynomial of degree at most degree; after rounding, both hold to rounding.

At n = 10 and degree 3 the end-point weights, at ahead 0, are (-56, 49, 64, 24, -36, -81, -76, 14, 224, 589)/715, and the one-bar forecast weights, at ahead 1, are (-12, 8, 13, 8, -2, -12, -17, -12, 8, 48)/30, each rounded to the nearest float64. Degree 0 gives n copies of the float64 nearest 1/n for every ahead. PolyFitWeights(2, 1, 0) is [0 1], and n = degree+1 interpolates the window, so that at ahead 0 the weights are 0, …, 0, 1.

The weights reverse with ahead: the reverse of PolyFitWeights(n, degree, ahead) is PolyFitWeights(n, degree, -(n-1)-ahead), bit for bit, when float64(1-n) - ahead, computed in float64, equals -(n-1)-ahead exactly, as it does for a dyadic ahead of moderate size, so the centre, ahead -(n-1)/2, gives symmetric weights. Ahead 0 gives the end-point Savitzky-Golay weights and, for odd n, ahead -(n-1)/2 the centred ones.

The result is nil, before anything is allocated, when degree is outside 0 to 4, when n < degree+1 or when ahead is NaN or infinite. It is also nil when a weight rounds to an infinity, which happens only from degree 2 on, for a huge ahead: PolyFitWeights(3, 2, 1e154) is finite and PolyFitWeights(3, 2, 2e154) nil, while degrees 0 and 1 give finite weights for every finite ahead. A call costs O(n*degree) big-rational operations, on numbers whose size grows with n, the degree and the number of bits of ahead, and returns a new slice every time, so n must be allocatable.

RepeatedMedianSlope

func (s Series) RepeatedMedianSlope(n int) Series

RepeatedMedianSlope returns Siegel’s repeated-median slope of the window of the last n bars against bar number, the median over the window’s bars of the median of each bar’s pairwise slopes to the other bars, a robust substitute for the least-squares slope in units of s per bar.

The window at bar t holds lags 0 to n-1, the current bar included, as w[k] for k from 0 to n-1 with w[n-1] the current bar, and the abscissa of w[k] is k. For each bar k the n-1 pairwise slopes to the other bars are the divisions (w[m]-w[k])/float64(m-k) for every m other than k, and the inner median M[k] is their median. The result is the outer median, the median of the n inner medians, plus 0, so that a zero result is +0 rather than a negative zero. At both levels the median of an odd count is the middle value and the median of an even count is the mid-average (a+b)/2 of the two middle values a and b, one addition and one division rounded as written, with no fallback when a+b overflows; the inner count n-1 is even exactly when n is odd and the outer count n is even exactly when n is even. Every order statistic is fixed as a number, so the result is the same whatever selection algorithm is used, the final plus 0 settling the sign of a zero.

The method is Siegel’s (Robust Regression Using Repeated Medians, Biometrika 69:1, 1982), which Meyers trades as the repeated median velocity (Dennis Meyers, “The Repeated Median Velocity Strategy”, Parts 1 and 2, Stocks & Commodities V31:12 (November 2013) and V31:13 (December 2013)). Its breakdown point is Siegel’s 50%, reached as the window grows: while fewer than (n-1)/2 of a window’s points are outliers, the result stays within the range of the slopes between the other points, so isolated outliers that pull the least-squares slope leave it alone. Meyers’s printed formula has two errors, both corrected here: its numerator reads price(t) where price(t-j) is meant, the slope between bars t-j and t-i, as his repeat lines and his Figure 2 show; and its indices run from 1 to N, which read literally as lags 1 to N would drop the current bar, whereas his worked example and Siegel number the N points of a fit 1 to N, read here as lags 0 to N-1, the current bar included. Neither Meyers nor his example fixes the median of an even count, the example giving 1 under any rule; (a+b)/2 is the correctly rounded midpoint and keeps negation and reversal exact, which a low or high median or the interpolated a+0.5*(b-a) does not. Meyers’s threshold rules, long above vup and short below -vdn, are a recipe over comparisons and no part of the method. Meyers’s 15-point example, the values 1 to 15 with 10 and 14 replaced by 18, gives exactly 1 at the last bar, the slope of the line through the other points, where Slope(15) gives 8/7 (about 1.1429).

The warm-up reads 0, not NaN: bars 0 to n-2 are 0 and the first valid index is n-1. The minimum length is n, and n = len(s) gives one valid bar, the last. For n < 2 or n > len(s), math.MaxInt and math.MinInt included, the result is len(s) zeros and only the result is allocated.

A valid bar is NaN while its window holds a NaN or an infinity, so a NaN or an infinity at index p gives NaN on [max(p, n-1), p+n-1] and nothing beyond, and later bars recover. On finite data the result is finite as long as every value of the window is below 4e307 in magnitude, where no slope and no mid-average overflows; beyond that a slope or a mid-average can overflow, the result can be an infinity, and a bar where an inner median is NaN, its two middle slopes being infinities of both signs, is NaN, which on finite data needs an odd n.

Each valid bar depends only on its own n values, so the result is window-local and start-invariant, bit for bit: dropping bars from the front of s leaves every later valid bar unchanged, and s[:j] gives the first j values. The method reads s only, so concurrent calls on a shared input are safe. The cost is O(n^2) expected time per valid bar, and O(n^2 log n) at worst, by selection rather than a full sort of each row of slopes; a call makes at most two allocations, the result and one scratch buffer of 2n-1 values, whatever len(s), and none per bar; there is no cap on n. Siegel’s repeated-median intercept and the Theil-Sen slope are not provided.

Bit for bit, on finite data below 4e307 in magnitude, n = 2 equals Delta(1) (bar 0 is 0 in both) apart from a -0 change of Delta, which reads +0 here. On every input, negating the series negates the result, as does reversing the window in time, a zero result being +0 either way; scaling by a power of two scales the result exactly when nothing overflows, underflows or is subnormal; and on a line whose values and slope are exact binary fractions, such as 1234.5 + 0.25 times t, every valid bar is the slope exactly, as is Slope(n) when the slope has few significant bits, as 0.25 has. To rounding, adding a constant leaves the result unchanged, adding c times the bar number adds c, scaling by any factor scales it, and RepeatedMedianSlope(3) equals Slope(3), the least-squares slope it robustly replaces.

Verification: two clean-room references, in Python and Julia, evaluate the same binary64 program with an explicit (a+b)/2 and agree with it bit for bit; the slope of scipy’s siegelslopes(y, arange(n)) is a cross-check on finite windows within the domain and free of subnormal values (scipy forms each slope in the opposite orientation, which gives the same bits once 0 is added to each slope, as here, and beyond the domain returns a finite value where this method overflows); R’s median is not (a+b)/2, so an R translation must code the mid-average by hand.

SlopePercent

func (s Series) SlopePercent(n int) Series

SlopePercent returns the n-bar least-squares slope as a percentage of the first price in its window, a scale-free trend rate per bar.

It is the percentage linear-regression slope of Markos Katsanos, Intermarket Trading Strategies (Wiley, 2008), the leg of his linear-regression-slope divergence (section 9.4, pages 127-128), in the form of his published code (Appendix A.1, page 298), “LinRegSlope(C,D1)/Abs(Ref(C,-D1+1))*100”. For j >= n-1 the output is signedPctOf(s.Slope(n)[j], s[j-n+1]): the ordinary least-squares slope of the window against bar number (Slope, in price units per bar), divided by the absolute value of the window’s first bar and then multiplied by 100, in his operation order. On a positive first bar it is (s.Slope(n)[j]/s[j-n+1])*100; on a negative first bar it divides by the absolute value, as his Abs does, so a rising negative series reads as rising. The unit is percent of the first price per bar.

The first valid index is n - 1 and the minimum length is n; every output before the first valid index is 0. An empty input, n < 2 or n > len(s) (math.MaxInt and math.MinInt included) returns a fresh all-zero series of len(s) in O(len(s)) time, without computing any slope.

It is window-local and start-invariant: each output depends only on the n values of its own window, so dropping bars from the front of s leaves every later output unchanged. Valid data: price levels, or any series whose windows start on a finite non-zero bar.

A NaN or an infinity at index p gives NaN on [max(p, n-1), p+n-1], and later outputs recover; a zero at p gives NaN at p+n-1 only, the end of the window that starts there, whose base is zero. In general a zero or non-finite first bar, a window holding a NaN or an infinity, or a non-finite result gives NaN. A flat window on a non-zero level gives exactly 0.

It is a portable trend filter, in percent per bar. Katsanos’s DAX system counts a 20-bar value below 0.2 percent per bar, with a congestion condition, as a trading range (“LRS20< .2”), and his FTSE system takes pullback buys only if the 80-bar value 15 bars earlier exceeds 0.03, which his code computes as “REF(100*LinRegSlope(C,80)/Ref(C,-79),-15) > .03”, equal to it in exact arithmetic on positive prices. Rafter’s moving slope rate of change (Stocks & Commodities 23:9), the (n+1)-bar slope divided by the price n bars back, is SlopePercent(n+1)/100 on positive prices: his quotient times 100 equals SlopePercent(n+1) bit for bit, but SlopePercent(n+1) divided by 100 recovers his quotient only in exact arithmetic. Katsanos’s soybean system normalises by the fitted line one bar before the window, not by the first price, so it is a different quantity.

Deviation from the MetaStock code: every output before the first valid index is 0, where MetaStock has no value.

The warm-up 0 fires Katsanos’s “LRS20< .2” trading-range rule on every bar before the first valid index, so trim (Window(-(n-1))) or mask (SetN(n-1, math.NaN())) the warm-up first.

Verification: two independent references, a Python reference in exact rational arithmetic and a Julia reference that also translates Katsanos’s “LinRegSlope(C,D1)/Abs(Ref(C,-D1+1))*100”, agree with it; it equals s.Slope(n)[i]/s[i-n+1]*100 bit for bit wherever the window’s first bar is positive and the result is finite; and the negation and power-of-two scaling identities hold bit for bit, apart from the sign of a zero and of a NaN: s.Mul(-1).SlopePercent(n) equals s.SlopePercent(n).Mul(-1), and a power-of-two scale of s that neither overflows nor underflows leaves the result unchanged.

Compositions

Chande’s forecast oscillator reads the regression line’s one-bar-ahead forecast, c.PolyFit(n, 1, 1), which equals tsf, the c.LinRegValue(n).Add(c.Slope(n)) of Meyers’s regression forecast on t6-regression-forecast.md below (TA-Lib’s TSF), to rounding. Each section gives its first valid index; prices are positive and finite and lengths are in bars.

Regression forecast oscillator %F (Chande)

Tushar Chande (“Forecasting Tomorrow’s Trading Day”, Stocks & Commodities V10:5, May 1992, pp. 220-224) fits “a simple linear regression” to five days and forecasts “the next day’s close with the slope m and constant c”, then defines “%F = ((Y - YForecast) / Y)100 where Y is the close for Day 6 and Y(Forecast) is the forecast for Day 6 from step 2 (from Day 5)”. He reads it with the fit’s r^2: “Use the forecasts only if r^2 is greater than 0.1”, and “A trend is in place if r^2 is greater than 0.6”. t6ForecastOscillator uses PolyFit(n, 1, 1), the least-squares line of the last n values evaluated one bar ahead (the editor’s sidebar builds the same cells in Excel with LINEST).

f := c.PolyFit(n, 1, 1).Shift(1).SetN(n, math.NaN())
pctF := c.Sub(f).Div(c).Mul(100)
r2 := c.R2(n).SetN(n-1, math.NaN())
fHigh := h.PolyFit(n, 1, 1).SetN(n-1, math.NaN())
fLow := l.PolyFit(n, 1, 1).SetN(n-1, math.NaN())

Inputs: the high h, low l and close c and the length n (5); pctF is in percent of the close, r2 runs from 0 to 1, and fHigh and fLow are the forecasts for bar t+1 made at bar t, Chande’s stop levels.
First valid index: n for pctF and n - 1 for r2, fHigh and fLow (5 and 4).
Anchors: on the sidebar’s OEX sheet (2 to 16 October 1991) every printed cell of rows 8 to 14 lies within 0.005 of the composition: the constant c.Average(5) less 3 times c.Slope(5), the slope, r^2 (0.72, 0.71, 0.00, 0.18, 0.56, 0.96, 0.95), the forecast c.PolyFit(5, 1, 1).Shift(1) and %F (0.01, 1.31, 0.57, 1.41, 0.99, -0.28); six cells are exact half-way ties at two decimals, which the sheet printed both ways. The %F numerator is a six-tap filter: c.Sub(c.PolyFit(5, 1, 1).Shift(1)) equals c.FIRSum([]float64{0.4, 0.1, -0.2, -0.5, -0.8, 1}) from bar 5 bit for bit (up to the sign of a zero). c.PolyFit(n, 1, 1) equals c.LinRegValue(n).Add(c.Slope(n)) to rounding. On 100 + 0.25t + q t^2 the miss c - f is q(n + 1)(n + 2)/6 at every valid bar, to rounding (0.07 at n = 5 and q = 0.01), so %F takes the sign of the acceleration, as in “The trend will usually change in the direction of %F”.
Traps: unmasked, the shifted forecast is 0 on bars 0 to n - 1 and %F reads 100 there; c.R2(n) and c.Slope(n) read 0 in their warm-up, which r2.LT(0.1) would take as no trend. A flat window gives r^2 and slope 0, and a zero close gives %F 0 through Div, where Excel divides by 0. The table’s bands are half-open (0.1<=r^2<0.3 is weak) while the rules say “greater than”: the gates are r2.GT(0.1) and r2.GT(0.6), the grade r2.GTE(0.1).Add(r2.GTE(0.3)).Add(r2.GTE(0.6)).
Variant: pandas-ta’s cfo is unlagged, today’s close against the forecast for tomorrow: c.Sub(c.PolyFit(n, 1, 1)).Div(c).Mul(100).
Missing piece: the stops “a few ticks above” fHigh and below fLow are orders: tier B (algolang).
Evidence: Chande’s charts of the OEX from October 1991 to January 1992 and of cash wheat in 1989; the sample plan’s data for 17 January 1992 is announced but not printed. Evidence, not a reference.
Pinned by TestRecipeT6ForecastOscillator; runnable example ExampleSeries_PolyFit_forecastOscillator.

The regression line’s one-bar forecast, the significance of its slope and its slope per average true range. Values are NaN before their first valid index, where the fit’s warm-up reads 0 or an average passes prices through; flags and states read 0 there, a state’s 0 being undetermined. Prices are positive and finite and lengths are in bars; reversal thresholds are fractions written as literals (0.01425 for 1.425%), rate-of-change thresholds percent; series are paired by position, so the caller aligns them by timestamp first.

Linear regression forecast with percentage reversal (Meyers)

Dennis Meyers (“Surfing The Linear Regression Curve With Bond Futures”, Stocks & Commodities V16:5, May 1998, pp. 209-218) follows “the plotted curve of pf”, the T-day least-squares line’s forecast of the next close, “Pf = a0 + b0*(T+1)”:“If pf has moved up by more than the percentage amount of pctup from the lowest low recorded in pf while short, then buy”, the mirror rule sells, and “if pf is greater than pf two days ago by the percentage amount %jmpup, we will go long”, with its mirror. t6RegressionForecast takes the forecast as LinRegValue plus Slope, TA-Lib’s TSF, runs his reversal with DirectionalChange and adds the jumps by running the difference machine on the log of the curve, where a jump bar’s threshold is 0 and so any move against the run reverses it.

tsf := c.LinRegValue(n).Add(c.Slope(n)).SetN(n-1, math.NaN())
state := tsf.DirectionalChange(up, down)
roc := tsf.RateOfChange(2)
hu := roc.GT(jumpUp).Mul(-1.0).Add(1.0).Mul(up).Log1p()
hd := roc.LT(-jumpDn).Mul(-1.0).Add(1.0).Mul(-down).Log1p().Mul(-1.0)
system := tsf.Log().DirectionalChangePointsScaled(hu, hd)

Inputs: positive closes; n (30; 33 in his first windows); up and down as fraction literals (his last window’s 0.01425 and 0.01425); jumpUp and jumpDn in percent (0.35 and 0.50).
First valid index: n - 1 for tsf, NaN before; state and system have no warm-up of their own and read 0 until their first confirmation (bars 65 and 57 on the synthetic path at n = 30): the masked bars are missing bars for both machines. They are the decisions after the close; Meyers fills at the next open.
Anchors: on Merrill’s sidebar data (1, 3, 2, 4, 3; “Fitting A Trendline By Least Squares”, in the same article) the slope is 0.5, the line at the fifth point 3.6 and the forecast 4.1, bit for bit (PolyFit(5, 1, 1) reads 4.1000000000000005); on Chande’s eleven closes the forecast equals c.PolyFit(5, 1, 1) to rounding (the %F section of this category pins his sheet); state equals a loop of the rule with the extreme restarted at the signal bar, and system the same loop with zero thresholds on jump bars; a jump exactly at its level does not fire.
Traps: the reduced forecast form printed on his p. 212 is a misprint: its weights, 2(3t - 2T - 4)/(T(T + 1)), sum to -(T + 5)/(T + 1) and give -10/3 on Merrill’s data, where the forecast is 41/10; the correct weights are 2(3t - T - 2)/(T(T - 1)). The Average form c.Average(n).Add(c.Slope(n).Mul(float64(n+1) / 2)), the PMA form and PolyFit equal the forecast only to 4 ulp, and the first two pass prices through in the warm-up, which the reversal then tracks. A run-time 1.35/100 is 0.013500000000000002, one ulp above 0.0135 (his third window’s pctdn), so write the fractions as literals. At a threshold tie the log differences of system and his fraction tests can split; on the synthetic paths every move stays at least 5.9e-5 from its threshold.
Variant: his earlier UD% rule (“UD% Price Change Indicator”, Stocks & Commodities V13, July 1995, pp. 298-303), the same reversal on the close, c.DirectionalChange(u, d).
Missing piece: an exact jump override needs per-bar fraction thresholds, DirectionalChangeScaled(theta Series, thetaDown ...Series), which the shipped DirectionalChangePointsScaled gives only in points.
Evidence: his walk-forward test (back-adjusted US Treasury bond futures, 1988 to February 1998, parameters re-optimised on five-year windows): 26 out-of-sample trades, 58% profitable and 55,875 dollars before costs, with in-sample trade lists that do not reconcile with his summaries. Evidence, not a reference.
Pinned by TestRecipeT6RegressionForecast; runnable example ExampleSeries_LinRegValue_meyersForecast.

Normalised slope: the t-statistic and R-squared significance (Siroky; Star; Karczewski)

Mike Siroky (“Normalized Slope As A Measure Of Trend”, Stocks & Commodities V38:7, June 2020, pp. 26-29) divides the slope by its standard error and notes that “this normalized slope is the same as the t-statistic for the slope”. Barbara Star (“Confirming Price Trend”, Stocks & Commodities V25:13, December 2007, pp. 72-78) reads R-squared against its 95% level: “a 10-period r-squared reaches a positive correlation at the 0.40 level, while a 50-period only needs to exceed a 0.08 level”, and 0.20 at 20. Jack Karczewski (“Identifying Market Trends”, Stocks & Commodities V16:4, April 1998, pp. 173-179) calls readings “of greater than 0.70” over 30 days a good fit. t6NormalisedSlope forms t as r sqrt((n - 2)/(1 - r^2)) from CorrelationTrend with Pow(-0.5), and gates its sign with R-squared.

r := c.CorrelationTrend(n)
t := r.Mul(r.Squared().Mul(-1.0).Add(1.0).Div(n-2).Pow(-0.5)).SetN(n-1, math.NaN())
r2 := c.R2(n).SetN(n-1, math.NaN())
trend := t.GT(0.0).Sub(t.LT(0.0)).Mul(r2.GT(crit))

Inputs: closes, n >= 3 (Siroky 26, Star 20, Karczewski 30) and crit: Star’s tc^2/(tc^2 + n - 2), with tc Student’s two-sided 5% value at n - 2 degrees of freedom, as a literal (0.3993 at n = 10, 0.1969 at 20, 0.0777 at 50), or Karczewski’s 0.7.
First valid index: n - 1 for t and r2, NaN before and on windows that hold a NaN; trend reads 0 there, and -0 where t is negative and R-squared at or below crit.
Anchors: Karczewski’s two windows, read off his Figures 1 and 2, give R-squared values that print as his 0.00014 and 0.95 (36/249023 and 1834854/1926557 exactly) and t values of 0.0636 and 23.669; t^2/(t^2 + n - 2) is c.R2(n) to 1e-15 and t has the sign of Slope on every valid bar; a perfect fit gives t = +Inf or -Inf, so the gates fire, and a flat window 0; only n - 2 degrees of freedom gives Star’s 0.40 at n = 10.
Traps: the dossier’s r.Mul(math.Sqrt(float64(n-2))).Div(r.Squared().Mul(-1.0).Add(1.0).Sqrt()) gives 0 at a perfect fit (Div’s zero rule), reading a perfect trend as none; Siroky counts n - 1 degrees of freedom (his 2.13 at n = 16 is Student’s t at 15), wrong for a slope; on price levels the textbook 5% level rejects on 55% to 82% of random-walk windows (n 10 to 50), so the levels are conventions, not p-values; Star’s eSignal slope code numbers x from 1 in one loop and from 0 in the other, which scales the slope by (n^2 - 1)/(n^2 + 11), and runs x back in time (0 is the latest bar), which her slope sidebar undoes by negating the result.
Variant: Karczewski’s runs, trending.SetN(n-1, math.NaN()).StateAge().Mul(trending) with trending := r2.GT(0.7), and his percentage error, c.Sub(fit).Div(fit).Mul(100.0) with fit := c.LinRegValue(n).
Variant: t without the ramp, sl.GT(0.0).Sub(sl.LT(0.0)).Mul(q.Mul(float64(n - 2)).Mul(q.Mul(-1.0).Add(1.0).Pow(-1.0)).Sqrt()) with sl := c.Slope(n) and q := c.R2(n), equal to rounding; for plotting, t.MaskNonFinite() turns the infinities into NaN.
Evidence: Katsanos’s R-squared system (“Which Trend Indicator Wins?”, Stocks & Commodities V34:11, October 2016, pp. 12-18) squares a 3-bar average of r, c.CorrelationTrend(18).Average(3).Squared(), not R-squared, and made 17 out-of-sample trades with a profit factor of 4.18 on one ETF, with parameters fitted in sample. Evidence, not a reference.
Pinned by TestRecipeT6NormalisedSlope; runnable example ExampleSeries_CorrelationTrend_normalisedSlope.

Linear trend per ATR (Masters; Aronson and Masters)

David Aronson and Timothy Masters’s manual for TSSB (Statistically Sound Machine Learning for Algorithmic Trading of Financial Instruments, 2014, p. 127) defines LINEAR PER ATR as “the slope of the line divided by the ATR”, fitted to “the log of the mean of the open, high, low, and close”, with “a small amount of compression and rescaling” that it does not print. Masters’s later public code (Statistically Sound Indicators for Financial Market Prediction, 2019, COMP_VAR.CPP, LINEAR TREND) fits the log close, divides by the mean log true range, multiplies by the fit’s R-squared (“Degrade the indicator if it is a poor fit”) and compresses with 100 Phi(x) - 50. t6LinearTrendATR is that code.

lc := c.Log()
atr := c.ATR(h.Log(), l.Log(), lc, m)
x := lc.Slope(n).Div(atr).Mul(scale).Mul(lc.R2(n)).SetN(max(n-1, m), math.NaN())
v := x.Mul(0.7071067811865476).Erf().Mul(50.0)

Inputs: positive highs, lows and closes; n and m (10 and 252 in his LINTRND_10); and scale, his Legendre normalisation, computed by the caller as math.Sqrt(float64(n*(n+1)) / float64(3*(n-1))) for n >= 3 (2.018433569398328 at n = 10, 3.0046260628866577 at 25) and the constant 1 / math.Sqrt2, 0.7071067811865476, at n = 2.
First valid index: max(n - 1, m), his front_bad, NaN before, where he writes 0.
Anchors: x and v equal a translation of his loop (the unit Legendre coefficient times 2 over n - 1 times his mean log true range, the R-squared and the exact normal distribution) to 1e-11; c.ATR(h.Log(), l.Log(), c.Log(), m) is his mean of log(max(H/L, H/C’, C’/L)), C’ the previous close, to rounding from bar m: a simple mean, not ATRWilder.
Traps: a run-time 1/math.Sqrt(2) is 0.7071067811865475, one ulp below the constant; where the average true range is 0 and the slope is not (n > m + 1, after a flat run), x and v read 0 by Div’s zero rule, where his 1e-60 guard gives +50 or -50; his own normal distribution is an approximation within 7.5e-8 (7.5e-6 on this scale); the TSSB manual’s form (the log of the bar’s mean price, no R-squared factor, an unprinted compression) is another indicator, and the “Squash transforms” of the Math introduction belong to TSSB’s separate SCALE and NORMALIZE options.
Variant: TSSB’s bare velocity, p.Slope(n).Div(p.ATR(h.Log(), l.Log(), c.Log(), m)) with p := o.SAverage(h, l, c).Log(), uncompressed.
Evidence: in the manual’s five-model subsampling check of its stand-alone system, whose second most used input is LIN_ATR_5, the out-of-sample profit factors are 1.069 to 1.088, with no cost model and no evidence for the indicator alone. Evidence, not a reference.
Pinned by TestRecipeT6LinearTrendATR.

Clifton Mitchell’s buy setup on the regression line of the lows, with an order that rests at the prior close until a later bar reaches it. Recipe values and levels are NaN before their first valid index (SetN(F, math.NaN())), a fill price is meaningful on fill bars only, and the setup, fill and book flags read 0 there: the regression line’s +0 warm-up is masked before a positive low is compared with it. Prices are positive and finite and the bars daily; series are paired by position, so the caller aligns them by timestamp first. The recipe reads the open o for its fill price only. A resting buy stop is touched when the high reaches it (GTE), as Mitchell’s “equal to or higher than the previous close” says, and fills at max(o, level), written o.Maximum(level): the level, or the open when the bar opens beyond it; the exits rest at levels known at the open and fill at o.Minimum(level). A high or low exactly at a level built with a decimal multiple of the ATR (the 1.2 of exit (c), the 0.25 of Acme V) is decided by binary rounding; the goldens run on prices in sixty-fourths. The backtests quoted are over data we do not have: evidence, not a reference. Counts on probe paths that the tests do not reproduce are illustrations, not pins.

Regression-of-lows setup with a pending close stop (Mitchell)

Clifton Mitchell (“Position Trading The S&P”, Stocks & Commodities V19:10, October 2001, pp. 26-33) sets the regression’s target “to predict for the current bar” and writes the rule as “If today’s low > yesterday’s linear regression of last X lows (buy setup), then buy at subsequent close stop”: “on any day following the buy setup, if the price was equal to or higher than the previous close, a buy signal would be triggered”, and the order “remained thus until the price reached the previous day’s close”. His second exit sells “if today’s low dropped below yesterday’s linear regression price (target), calculated from the last X days’ lows”, over three days. t6RegressionLows returns the setup, the fill of the pending order (kept live by ValueWhen over a code of setup and trigger bars), its price px and the book pos with that exit.

lrPrev := l.LinRegValue(n).Shift(1).SetN(n, math.NaN())
setup := l.GT(lrPrev)
trig := h.GTE(c.Shift(1))
code := setup.Mul(2).Add(trig)
armed := code.ValueWhen(code).Shift(1).GTE(2)
fill := armed.Mul(trig)
px := o.Maximum(c.Shift(1))
lr3 := l.LinRegValue(3).Shift(1).SetN(3, math.NaN())
pos := c.Position(fill, l.LT(lr3))

Inputs: daily opens, highs, lows and closes; n, the entry regression length: 9 with exit (a) (“the entry linear regression lookback period optimized at nine days” for the low stop), 6 with exit (b), whose own length is 3, and 6 with exit (c), an 8-day ATR times 1.2 subtracted from the highest high of 4 days. His market is the S&P 500 index.
First valid index: n for setup (bars 0 to n - 1 read 0 through the NaN mask) and n + 1 for fill; px is meaningful on fill bars only (on bar 0 it is max(o[0], +0), the open); pos reads 0 until the first fill. The level of exit (c) is NaN before bar 8 (ATR(8) passes the true range through on bars 0 to 6). trig reads 1 on bar 0 against the +0 pad, so bar 0 is an event (code 1) and code.ValueWhen(code) has no NaN on any bar; armed is unaffected: it reads the Shift pad on bar 0, and a trigger-only bar arms nothing. The SetN(3, math.NaN()) of lr3 is a no-op on positive prices (the line’s +0 warm-up, shifted, is below every positive low): the test asserts the masked and unmasked forms equal.
Anchors: code is 2 or 3 on setup bars, 1 on trigger-only bars and 0 elsewhere, and armed is 1 when the latest bar with a code, as of the bar before, was a setup: the order lives from the bar after a setup until the first bar whose high reaches the prior close, and a setup bar that also triggers counts as a setup, the order being placed after its close. fill equals a loop of the pending order bit for bit (on the 320-bar group and three 3,000-bar walks, 151 and 1,452 to 1,504 fills). The order is a function of prices alone: a trigger bar while long consumes it, and Position ignores a fill while long and tests no exit on the entry bar (the source is silent on orders armed during a trade). l.LinRegValue(n) equals l.Average(n).Add(l.Slope(n).Mul(float64(n-1) / 2)) bit for bit from bar n - 1 (+0 before, where Average passes the lows through), and l.LinRegValue(3) equals the dossier’s l.SMA(3).Add(l.Slope(3)) from bar 2. Highs exactly at the prior close and a low equal to lr3 that holds are planted on the fixture; the strict setup is pinned by a constructed case. The trade-state page t2-latched-entries.md (“Stop entry above the prior close (Connors)”) keeps its order for one bar, setup.Shift(1).Mul(h.GTE(lvl)), the special case of this one that cancels after a bar when the level is c.Shift(1); Meyers’s forecast on this category’s t6-regression-forecast.md is the same line one bar ahead, LinRegValue(n).Add(Slope(n)).
Traps: l.PMA(n) is not the end point but half a bar past it (PolyFit(n, 1, 0.5) to rounding, and it passes the lows through on bars 0 to n - 2): the setup built on it differs on about 1 bar in 15 (n 6) and 1 in 20 (n 9) of 20,000-bar probe walks. Unmasked, l.GT(l.LinRegValue(n).Shift(1)) fires on every warm-up bar of a positive series. The BarsSince form setup.BarsSince().Shift(1).LTE(trig.BarsSince().Shift(1)) is not the order: it reads 0 until the first trigger and compares two pads on bar 0. The dossier’s claims that the pending order needs an event layer and that nseries lacks a named regression end point are stale (ValueWhen keeps the order, and LinRegValue is the end point), and its close test for exit (c) with today in both windows, c.LT(h.Highest(4).Sub(c.ATR(h, l, c, 8).Mul(1.2))), is a materially different rule: fewer trades and more time in the market than the stop on probe walks (about 2,400 trades and 56% to 59% of bars, against 3,400 to 3,600 trades and 44% to 45%; not pinned).
Variant: Acme V (Mark R. Conway and Aaron N. Behle, Professional Stock Trading: System Design and Automation, Acme Trader, 2002, ch. 6, pp. 123-134, Example 6.1): with atr := c.ATR(h, l, c, 20), the setup l.GT(l.PolyFit(5, 1, 1)) (the code’s LinearRegValue(Low, 5, -1): the five lows ending today, projected one bar ahead and not shifted) times h.LT(h.Lowest(4).Shift(1)), h.Sub(l).LTE(atr.Mul(rf)) (RangeFactor rf 1.0), l.GT(c.Average(50)), h.GT(l.Shift(1)) and the caller’s 0/1 weekday mask wd, all masked with SetN(49, 0), then a one-bar buy stop at c.Add(atr.Mul(0.25)).Shift(1); its departures in one line: ATR(20) for the unpublished Volatility(20), the caller’s Monday-or-Tuesday mask on the setup bar (fills on Tuesday or Wednesday), the code’s <= range test where the rules box prints <, and ch. 6’s level, the close plus 0.25 ATR, where ch. 11 prints Close + EntryFactor. The prose’s window (four lows ending yesterday, l.PolyFit(4, 1, 1).Shift(1)) and the dossier’s (l.PolyFit(5, 1, 1).Shift(1)) are not the code’s, and the dossier’s differs from it on 32% of probe bars; PolyFit(5, 1, 1) is Meyers’s forecast quantity to rounding (at most 3 ulp).
Members:
- exit (a), the low stop, with n 9: c.Position(fill, l.LTE(l.Shift(1))), filled at o.Minimum(l.Shift(1))
- exit (c), with n 6: lvl := h.Highest(4).Sub(c.ATR(h, l, c, 8).Mul(1.2)).Shift(1).SetN(8, math.NaN()), l.LT(lvl)
Evidence: his Figure 2 (the S&P 500 index from January 1995 to May 2001, 0.5 point a trade for slippage and commissions, the lengths optimised on the same data, in sample) gives 341, 345 and 224 trades for exits (a), (b) and (c), profit factors of 2.03, 2.07 and 2.22 and 1201.47, 1217.61 and 1206.67 points, against 797 for buy and hold, and he warns that “the profits may be more of a result of the overall market rise”; Conway and Behle’s Tables 6.1 and 6.2 (over a thousand stocks, March 2000 to March 2002, costs not stated) print profit factors of 1.75 and 1.20. Evidence, not a reference.
Pinned by TestRecipeT6RegressionLows; runnable example ExampleSeries_LinRegValue_regressionOfLows.

Fractional Differencing

FracDiff

func (s Series) FracDiff(d float64, n int) Series

FracDiff returns the fixed-width fractional difference of order d of s, the truncated binomial expansion of the operator (1-B)^d, B the backshift operator, applied to s over a window of n weights, so that out[t] is the sum of w_k s[t-k] for k = 0 to n-1. It is the fixed-width window form of López de Prado, Advances in Financial Machine Learning (2018), here “AFML”, section 5.5.2 and Snippet 5.3, on the operator of Hosking (1981) and of Peters, Fractal Market Analysis (1994), chapter 13.

The result is s.FIRSum(FracDiffWeights(d, n)) bit for bit: each product is rounded on its own, the terms are summed oldest first from +0 and 0 is added at the end, so a zero result is +0 and every platform gives the same bits.

n is the number of weights, covering lags 0 to n-1: it is the window of the Python package fracdiff and AFML’s “width + 1”, since Snippet 5.3 counts lags. AFML’s and mlfinpy’s threshold form is s.FracDiff(d, FracDiffWidth(d, tau)), which is AFML’s fracDiff_FFD to rounding on gap-free finite data; FracDiffWidth’s 0 is an invalid n here, so a failure of the threshold rule shows as len(s) zeros.

Any finite d is accepted. d = 0 gives the input plus 0 from bar n-1, the other weights being -0. d = 1 is the first difference, Delta(1): with n = 2 it matches over the whole series except for the sign of a zero, and with n > 2 it matches from bar n-1 on finite data. d = 2 with n = 3 is the second difference, s.Delta(1).Delta(1) from bar 2, to rounding. A negative d is truncated fractional integration: d = -1 gives the n-bar window sum, Sum(n) from bar n-1 to rounding, and d = -2 the weights 1, 2, 3 and so on.

The method is linear in s and takes no logarithm: pass log prices yourself, as in Log().FracDiff(d, n), for example c.Log().FracDiff(d, n), which is what AFML and Peters intend, while Kaabar, Deep Learning for Finance (2024), applies fracdiff to price levels.

A constant input c gives cw_{n-1}(d-1), that is c times FracDiffWeights(d-1, n)[0], to rounding, and that is 0 only for an integer d from 1 to n-1, so the output on log prices is “level plus noise” and not a zero-mean return. There is no flat-window rule: a flat window of x gives xw_{n-1}(d-1) to rounding.

The warm-up reads 0, not NaN: bars 0 to n-2 are 0 and the first valid index is n-1. Trim it with Window(-(n-1)) for n >= 2, since Window(0) is empty, or mask it with SetN(n-1, math.NaN()). n = 1 gives each value plus 0. n < 1, n > len(s), a non-finite d and overflowing weights give len(s) zeros, and n > len(s) is checked before any weights are built, so a huge n allocates only the zeros.

Each valid bar costs O(n) and depends only on its own window, so a prefix of s gives a prefix of the result, and the weights cost O(n) memory besides the result, a new series of len(s) values; s is not modified.

AFML picks the smallest d whose output passes an ADF test (Snippet 5.4: statsmodels adfuller, one lag, a constant, tau = 0.01), and nseries has no ADF test with lagged differences (DickeyFuller is the form without lags). Peters’s d = H - 0.5 (his equation 13.10) is theory for the return process, so on log prices the implied order would be H + 0.5, and Hurst is not wired in.

FracDiffWeights

func FracDiffWeights(d float64, n int) []float64

FracDiffWeights returns the n weights of the fixed-width fractional difference of order d, oldest first: index n-1 holds w_0 = 1, the current bar’s weight, and index 0 holds w_{n-1}, the oldest bar’s weight, so that s.FIRSum(FracDiffWeights(d, n)) is s.FracDiff(d, n).

In closed form w_k = (-1)^k C(d, k), the product of (i - d)/(i + 1) for i = 0 to k-1. They are built instead by the recursion of Snippet 5.1 of López de Prado’s Advances in Financial Machine Learning (2018, “AFML” below and in FracDiff), with w_0 = 1 and, for k = 1 to n-1,

w[k] = -w[k-1] / float64(k) * (d - float64(k) + 1)

which divides before it multiplies (mlfinpy multiplies first) and forms the last factor as (d - float64(k)) + 1. Nothing is added to a product, so no compiler fuses the recursion and every platform gives the same weights. Because the last factor is rounded, w_1 is -d only to rounding: FracDiffWeights(0.1, 2) is [-0.09999999999999998 1], and a positive d up to 2^-54, about 5.6e-17, or a negative d down to -2^-53, about -1.1e-16, gives w_1 = -0. d = 0 gives 1 and then -0s, so FracDiffWeights(0, 4) is [-0 -0 -0 1]. A positive integer d gives the signed binomial row to rounding, exactly for d up to 28, and then exact zeros: FracDiffWeights(1, 4) is [0 0 -1 1] with +0s and FracDiffWeights(2, 5) is [-0 -0 1 -2 1]. From d = 29 a quotient such as 406/3 rounds, so FracDiffWeights(29, 30) holds w_3 = -3654.0000000000005 and ends in w_29 = -1.0000000000000007, while the zeros after k = d stay exact. d = 0.5 gives the dyadic weights -Catalan(k-1)/2^(2k-1), exact through k = 31. d = -1 gives 1 to rounding, exactly 1 through k = 48 with w_49 = 0.9999999999999999, and d = -2 gives exactly 1, 2, 3 and so on, so FracDiffWeights(-2, 6) is [6 5 4 3 2 1]. In exact arithmetic the weights sum to w_{n-1}(d-1), the last weight of order d-1.

The result is nil when n < 1 and when d is NaN or an infinity, in both cases before anything is allocated, and nil when any w_k for k = 1 to n-1 is not finite, at the first such weight. Any finite d is otherwise valid, negative and large values included. Overflow happens only at large |d| or large n. A positive d above about 1029.33 overflows once n reaches the hump of the weights near k = d/2, so FracDiffWeights(1030, 500) is finite and FracDiffWeights(1030, 501) is nil, while the integer d = 1029 peaks at 1.43e308 and is finite at every n, reaching an exact zero at k = 1030. For d < -1 the weights grow like k^(-d-1), so d = -1030 overflows from n = 304 and d = -100 from n = 48832, while d = -3 would need n near 2e154.

The slice is newly allocated on every call, and a finite d with a huge n allocates n values, so n must be allocatable, as for the package’s other weight builders such as ALMAWeights and HannWeights.

FracDiffWidth

func FracDiffWidth(d, tau float64) int

FracDiffWidth returns the number of weights that the threshold rule of AFML (see FracDiff) keeps for the fractional order d and the threshold tau, counted as mlfinpy’s get_weights_ffd counts it: the first k >= 1 at which the absolute value of w_k falls below tau, so that the weights 0 to k-1 are kept and a weight equal to tau is kept. It is the n of s.FracDiff(d, FracDiffWidth(d, tau)), which is AFML’s fracDiff_FFD(series, d, thres=tau) to rounding on gap-free finite data. mlfinpy multiplies before it divides, so its count can differ from this one where a weight lies within rounding of tau.

The domain is 0 < tau < 1, since w_0 = 1 is always kept. The result is 0 for any other tau, including 0, -0, negative values, 1, values above 1, NaN and the infinities, for a NaN or infinite d, for a weight that is not finite, and for a threshold that no weight up to k = fracDiffMaxWidth, inclusive, falls below; the cap is 1 << 20, that is 1,048,576. A crossing at exactly k = fracDiffMaxWidth returns fracDiffMaxWidth and one at fracDiffMaxWidth + 1 returns 0. FracDiff reads 0 as an invalid n and so gives zeros, which makes the failure visible, while mlfinpy instead truncates the width to the series length.

A tau above the absolute value of w_1 gives 1, as does d = 0 and, since w_1 is -d only to rounding, a d whose absolute value is below tau by more than that rounding. Equality keeps, so FracDiffWidth(0.5, 0.0625) is 4 and FracDiffWidth(0.5, 0.5) is 2. FracDiffWidth(0.4, 1e-5) is 1458, the window behind AFML’s Figure 5.4 (d = 0.4, tau = 1e-5), a count AFML does not print, and the largest count at tau = 1e-5 over d >= 0 is 4093, near d = 0.11, with 4076 at d = 0.1, 4078 at d = 0.12 and 929 at d = 0.01. For d in (0, 1) the largest counts per threshold are 12 at 1e-2, 74 at 1e-3, 528 at 1e-4, 4093 at 1e-5, 33,244 at 1e-6 and 279,095 at 1e-7, and for d >= 1 no count exceeds 1032 (an integer d gives d + 1, up to 1030 at d = 1029, and overflow ends the range just above), so the cap covers every tau >= 1e-7 for d >= 0, while tau = 1e-8 at d = 0.1 would need 2,174,855 weights and gives 0. An integer d from 0 to 1029 gives d + 1, and 1030 and above overflow and give 0. FracDiffWidth(0.9, 1e-3) is 12, mlfinpy’s own test. For d in (-1, 0) the weights decay like k^(-d-1), so d = -0.1 needs 29,411 at tau = 1e-5, while d = -0.5 would need about 3e9 and gives 0. For d <= -1 the result is 0 because the weights do not fall below tau, and d = -1’s weights stay within 1.5e-13 of 1 up to the cap, so only a tau that close to 1 is crossed: FracDiffWidth(-1, 0.9999999999999999) is 51 and FracDiffWidth(-1, 1 - 1e-14) is 15064. mlfinpy’s order keeps those weights at exactly 1, so it finds no crossing there.

FracDiffWidth allocates nothing and costs at most 1 << 20 steps.

Maths Moving averages and filters