Replacing magic numbers in Ehlers MAMA and HTIT

[mihakralj]

In the original Ehlers code for MAMA and HTIT, the magic numbers 0.0962 and 0.5769 are used to apply a Finite Impulse Response (FIR) filter to theoretical underlying IIR.

I hate poor rounded/truncated magic numbers, so I spent way too much time with Wolfram Alpha to find what actual ratios were used by Ehlers.

Passband Optimized Weights

• From: 0.0962 and 0.5769
• To: 5/52 and 15/26
• These coefficients define the Hilbert Transformer's frequency response. Using the exact fractions ensures the filter's gain is normalized and uses high precision. If you look at the relationship, 15/26 is exactly 6*(5/52). Using fractions prevents "drift" in the phase calculation caused by cumulative rounding errors in long time-series datasets. I know you cherish precision and hate drift as much as I do. :-)

Gain Regression Constants

The variables adj (adjustment) used 0.075 and 0.54 to correct the amplitude based on the detected period.

• From: 0.075 (Slope) and 0.54 (Intercept)
• To: 3/40 and 27/50
• These numbers come from a linear regression (y = mx + b) designed to compensate for the HT's inherent gain. 3/40 (0.075) is the specific rate at which the signal amplitude must be scaled relative to the period length. 27/50 (0.54) is the base offset.

Using ratios makes the origin of magic numbers in MAMA/HTIT clear: they are derived from geometric properties of the sine waves being analyzed, rather than being "tuned by hand" values.

Suggested change of your code in HTIT (similar would apply for MAMA):

namespace Skender.Stock.Indicators;

/// <summary>
/// Provides methods for calculating the Hilbert Transform Instantaneous Trendline (HTL) indicator.
/// </summary>
public static partial class HtTrendline
{
    /// <summary>
    /// Converts a list of time-series values to Hilbert Transform Instantaneous Trendline (HTL) results.
    /// </summary>
    /// <param name="source">The list of time-series values to transform.</param>
    /// <returns>A list of HTL results and smoothed price.</returns>
    /// <exception cref="ArgumentNullException">Thrown when the source list is null.</exception>
    public static IReadOnlyList<HtlResult> ToHtTrendline(
        this IReadOnlyList<IReusable> source)
    {
        // prefer HL2 when IQuote
        IReadOnlyList<IReusable> values
            = source.ToPreferredList(CandlePart.HL2);

        // initialize
        int length = values.Count;
        List<HtlResult> results = new(length);

        double pr = new double[length]; // price
        double sp = new double[length]; // smooth price
        double dt = new double[length]; // detrender
        double pd = new double[length]; // period

        double q1 = new double[length]; // quadrature
        double i1 = new double[length]; // in-phase

        double q2 = new double[length]; // adj. quadrature
        double i2 = new double[length]; // adj. in-phase

        double re = new double[length];
        double im = new double[length];

        double sd = new double[length]; // smooth period
        double it = new double[length]; // instantaneous trend (raw)

        // Fractional Ratios for Filter Coefficients and Linear Regressions
        const double CoeffA = 5d / 52d;    // Primary outer weights (0.09615...)
        const double CoeffB = 15d / 26d;   // Secondary inner weights (0.57692...)
        const double GainSlope = 3d / 40d; // Gain correction slope (0.075)
        const double GainInt = 27d / 50d;  // Gain correction intercept (0.54)

        // roll through source values
        for (int i = 0; i < length; i++)
        {
            IReusable s = values[i];
            pr[i] = s.Value;

            if (i > 5)
            {
                // Amplitude correction derived from previous period to maintain unity gain
                double adj = (GainSlope * pd[i - 1]) + GainInt;

                // 4-tap low-pass smoothing FIR filter
                sp[i] = ((4d * pr[i]) + (3d * pr[i - 1]) + (2d * pr[i - 2]) + pr[i - 3]) / 10d;

                // Improved Hilbert Transformer (7-tap FIR)
                dt[i] = ((CoeffA * sp[i]) + (CoeffB * sp[i - 2]) - (CoeffB * sp[i - 4]) - (CoeffA * sp[i - 6])) * adj;

                // Decompose into In-Phase and Quadrature components
                q1[i] = ((CoeffA * dt[i]) + (CoeffB * dt[i - 2]) - (CoeffB * dt[i - 4]) - (CoeffA * dt[i - 6])) * adj;
                i1[i] = dt[i - 3]; // 3-bar center-tap delay to align phases

                // Advance phases by 90 degrees to create analytic complex signal
                double jI = ((CoeffA * i1[i]) + (CoeffB * i1[i - 2]) - (CoeffB * i1[i - 4]) - (CoeffA * i1[i - 6])) * adj;
                double jQ = ((CoeffA * q1[i]) + (CoeffB * q1[i - 2]) - (CoeffB * q1[i - 4]) - (CoeffA * q1[i - 6])) * adj;

                // Phasor addition for temporal averaging
                i2[i] = i1[i] - jQ;
                q2[i] = q1[i] + jI;

                // Smoothing I and Q components (1/5 and 4/5 EMA)
                i2[i] = (0.2d * i2[i]) + (0.8d * i2[i - 1]);
                q2[i] = (0.2d * q2[i]) + (0.8d * q2[i - 1]);

                // Homodyne Discriminator: complex multiplication by previous conjugate
                re[i] = (i2[i] * i2[i - 1]) + (q2[i] * q2[i - 1]);
                im[i] = (i2[i] * q2[i - 1]) - (q2[i] * i2[i - 1]);

                re[i] = (0.2d * re[i]) + (0.8d * re[i - 1]);
                im[i] = (0.2d * im[i]) + (0.8d * im[i - 1]);

                // Calculate cycle period based on phase rotation (2*PI radians = 360 degrees)
                pd[i] = im[i]!= 0 && re[i]!= 0
                   ? 2d * Math.PI / Math.Atan(im[i] / re[i])
                    : 0d;

                // Constraint logic using exact ratio thresholds
                pd[i] = pd[i] > (3d / 2d) * pd[i - 1]? (3d / 2d) * pd[i - 1] : pd[i]; // 1.5x limit
                pd[i] = pd[i] < (67d / 100d) * pd[i - 1]? (67d / 100d) * pd[i - 1] : pd[i]; // 0.67x limit
                pd[i] = pd[i] < 6d? 6d : pd[i];
                pd[i] = pd[i] > 50d? 50d : pd[i];

                // Recursive smoothing of the period measurement
                pd[i] = (0.2d * pd[i]) + (0.8d * pd[i - 1]);
                sd[i] = (33d / 100d * pd[i]) + (67d / 100d * sd[i - 1]);

                // Calculate the Instantaneous Trendline over the Dominant Cycle
                int dcPeriods = (int)(double.IsNaN(sd[i])? 0 : sd[i] + 0.5d);
                double sumPr = 0;
                for (int d = i - dcPeriods + 1; d <= i; d++)
                {
                    if (d >= 0)
                    {
                        sumPr += pr[d];
                    }
                    else
                    {
                        dcPeriods--;
                    }
                }

                it[i] = dcPeriods > 0? sumPr / dcPeriods : pr[i];

                // Final indicator compilation
                results.Add(new(
                    Timestamp: s.Timestamp,
                    DcPeriods: dcPeriods > 0? dcPeriods : null,
                    Trendline: i >= 11 
                     ? (((4d * it[i]) + (3d * it[i - 1]) + (2d * it[i - 2]) + it[i - 3]) / 10d).NaN2Null()
                      : pr[i].NaN2Null(),
                    SmoothPrice: (((4d * pr[i]) + (3d * pr[i - 1]) + (2d * pr[i - 2]) + pr[i - 3]) / 10d).NaN2Null()));
            }
            else
            {
                // Initialization sequence
                results.Add(new(
                    Timestamp: s.Timestamp,
                    DcPeriods: null,
                    Trendline: pr[i].NaN2Null(),
                    SmoothPrice: null));

                pd[i] = 0;
                sp[i] = 0;
                dt[i] = 0;
                i1[i] = 0;
                q1[i] = 0;
                i2[i] = 0;
                q2[i] = 0;
                re[i] = 0;
                im[i] = 0;
                sd[i] = 0;
            }
        }

        return results;
    }
}

[DaveSkender]

Omg, I tried to figure this out a while ago and gave up!
If I recall, I concluded that these constants were an estimate based on a curve fitting value.
I'll take a look

[mihakralj]

The empirical coefficients 0.0962, 0.5769, 0.075, and 0.54 are discrete approximations of rational fractions derived from the optimization of a band-limited Hilbert Transformer and its associated gain-compensation model.

FIR Filter Weights: 5/52 and 15/26

The "Improved Hilbert Transformer" is a 7-tap anti-symmetric Finite Impulse Response (FIR) filter. Its objective is to implement a wideband 90-degree phase shifter for financial time series while providing rejection of high-frequency noise.

Mathematical Derivation of the 1:6 Ratio:

The transfer function in the Z-domain for an anti-symmetric FIR filter of this structure is:

H(z) = a + bz^-2 - bz^-4 - a*z^-6

John Ehlers optimized the response for market cycles typically observed between 10 and 40 bars. Through ripple minimization analysis, it was determined that the optimal weight relationship between the outer taps (a) and inner taps (b) follows a 1:6 ratio.

Dividing the standard implementation values confirms this relationship:

0.5769 / 0.0962 = 5.9968 (approximately 6.0)

Exact Fractional Rationalization:

To achieve high precision and eliminate rounding drift, these coefficients are expressed as rational fractions with a common denominator of 52:

• Coefficient A (Primary): 5/52 = 0.0961538... (rounded to 0.0962 in standard code).
• Coefficient B (Secondary): 30/52 = 15/26 = 0.5769230... (rounded to 0.5769 in standard code).

The use of "gapped" indices (using price data at i, i-2, i-4, and i-6) acts as an integrated filter that provides approximately -8 dB of rejection for 5-bar noise cycles, preventing jitter in the phase calculation.

Amplitude Correction: 3/40 and 27/50

A truncated Hilbert Transformer behaves like a numerical differentiator, meaning its output magnitude is not flat across the frequency spectrum. The amplitude exhibits a roll-off (attenuation) that is proportional to the cycle period.

Linear Regression Compensation Model:

To maintain unity gain, an amplitude correction factor is applied. This factor is a first-order linear approximation of the filter's gain error G(T) over the dominant cycle range (10 to 50 bars).

The linear model is defined as:

Correction(T) = m * T + b

Through a least-squares fit of the filter's theoretical response, the following coefficients were established:

• Slope (m): 0.075 = 75/1000 = 3/40.
• Intercept (b): 0.54 = 54/100 = 27/50.

This compensation ensures the In-Phase (I) and Quadrature (Q) components remain amplitude-balanced, which is a prerequisite for calculating the instantaneous phase: Phase = arctan(Q/I).

[DaveSkender]

@mihakralj this inspired me to pull out a couple of my old college books for a trip down memory lane:

• Applied Numerical Analysis has more relevant stuff for handling precision, but is word soup of _foundational computational methods. I was surprised to find a section on SIMD, apparently been around since the 1960's! Still reading ...
• Intro to Electrical Engineering has more on signal processing. I was looking at for more math around the wideband pass filters that Ehlers uses for his basis. Not finding much usefulness there.

Fun reads nonetheless.

[DaveSkender]

Ehlers' coefficients were an estimate based on a curve fitting value 👍

I did read more in Ehlers' Rocket Science for Traders book. I think you're over-imagining these fraction, but am still looking at some additional texts to find more of a mathematical basis.

is what he tweaked manually to get the truncated frequency band filter:

to look more flat like:

[DaveSkender]

This phase (waves) and Hilbert transform Wikipedia pages is a good supplemental on some of the elements of band pass filters.