using System.Runtime.CompilerServices; using System.Runtime.InteropServices; namespace QuanTAlib; /// /// PTA: Ehlers Precision Trend Analysis /// /// /// /// Dual highpass filter bandpass approach for near-zero-lag trend extraction. /// Applies two 2-pole Butterworth highpass filters with different cutoff periods /// to the same input, then subtracts: Trend = HP(longPeriod) - HP(shortPeriod). /// This preserves cyclic components between shortPeriod and longPeriod bars, /// producing a zero-centered trend indicator with near-zero lag. /// /// /// Algorithm based on: John F. Ehlers, "Precision Trend Analysis," TASC September 2024. /// /// /// Complexity: O(1) per bar — two IIR filter evaluations + subtraction. /// /// [SkipLocalsInit] public sealed class Pta : AbstractBase { // ── HP1 (long-period) coefficients ── private readonly double _c1L, _c2L, _c3L; // ── HP2 (short-period) coefficients ── private readonly double _c1S, _c2S, _c3S; [StructLayout(LayoutKind.Auto)] private record struct State { // HP1 (long-period highpass) public double Hp1; public double Hp1_1; // HP2 (short-period highpass) public double Hp2; public double Hp2_1; // Source history (shared by both filters) public double Src1; public double Src2; public int Count; public static State New() => new() { Hp1 = 0, Hp1_1 = 0, Hp2 = 0, Hp2_1 = 0, Src1 = 0, Src2 = 0, Count = 0 }; } private State _state; private State _p_state; /// Long-period cutoff for the first highpass filter. public int LongPeriod { get; } /// Short-period cutoff for the second highpass filter. public int ShortPeriod { get; } /// /// Initializes a new instance of the class. /// /// Long-period HP cutoff. Default is 250 (~1 year daily). /// Short-period HP cutoff. Default is 40 (~2 months daily). /// /// Thrown when longPeriod < 3, shortPeriod < 2, or longPeriod <= shortPeriod. /// public Pta(int longPeriod = 250, int shortPeriod = 40) { ArgumentOutOfRangeException.ThrowIfLessThan(longPeriod, 3, nameof(longPeriod)); ArgumentOutOfRangeException.ThrowIfLessThan(shortPeriod, 2, nameof(shortPeriod)); if (longPeriod <= shortPeriod) { throw new ArgumentOutOfRangeException(nameof(longPeriod), $"longPeriod ({longPeriod}) must be greater than shortPeriod ({shortPeriod})."); } LongPeriod = longPeriod; ShortPeriod = shortPeriod; // Precompute HP coefficients for long period ComputeHpCoefficients(longPeriod, out _c1L, out _c2L, out _c3L); // Precompute HP coefficients for short period ComputeHpCoefficients(shortPeriod, out _c1S, out _c2S, out _c3S); Name = $"PTA({longPeriod},{shortPeriod})"; WarmupPeriod = longPeriod; _state = State.New(); _p_state = _state; } /// /// Initializes a new instance of the class with a publisher source. /// public Pta(ITValuePublisher source, int longPeriod = 250, int shortPeriod = 40) : this(longPeriod, shortPeriod) { source.Pub += (object? _, in TValueEventArgs args) => Update(args.Value, args.IsNew); } /// /// Computes 2-pole Butterworth highpass filter coefficients. /// [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void ComputeHpCoefficients(int period, out double c1, out double c2, out double c3) { double f = 1.414 * Math.PI / period; double a1 = Math.Exp(-f); double b1 = 2.0 * a1 * Math.Cos(f); c2 = b1; c3 = -(a1 * a1); c1 = (1.0 + c2 - c3) * 0.25; } public override bool IsHot => _state.Count >= 2; /// Primes the indicator with historical data. public override void Prime(ReadOnlySpan source, TimeSpan? step = null) { foreach (double v in source) { Update(new TValue(DateTime.MinValue, v), isNew: true); } } [MethodImpl(MethodImplOptions.AggressiveInlining)] public override TValue Update(TValue input, bool isNew = true) { if (isNew) { _p_state = _state; } else { _state = _p_state; } double src = input.Value; ref State s = ref _state; double result; if (s.Count < 2) { // Bootstrap: not enough bars for HP differentiation if (s.Count == 0) { s.Src1 = src; s.Src2 = src; } else { s.Src2 = s.Src1; s.Src1 = src; } s.Count++; result = 0.0; } else { // 2nd-order difference: src - 2*src1 + src2 double diff = src - 2.0 * s.Src1 + s.Src2; // HP1 (long period): hp1 = c1L*diff + c2L*hp1 + c3L*hp1_1 double hp1 = Math.FusedMultiplyAdd(_c1L, diff, Math.FusedMultiplyAdd(_c2L, s.Hp1, _c3L * s.Hp1_1)); // HP2 (short period): hp2 = c1S*diff + c2S*hp2 + c3S*hp2_1 double hp2 = Math.FusedMultiplyAdd(_c1S, diff, Math.FusedMultiplyAdd(_c2S, s.Hp2, _c3S * s.Hp2_1)); // Trend = HP1 - HP2 (bandpass between shortPeriod and longPeriod) result = hp1 - hp2; // Update HP state s.Hp1_1 = s.Hp1; s.Hp1 = hp1; s.Hp2_1 = s.Hp2; s.Hp2 = hp2; // Update source history s.Src2 = s.Src1; s.Src1 = src; s.Count++; } Last = new TValue(input.Time, result); PubEvent(Last, isNew); return Last; } /// Updates with a full TSeries and returns results. [MethodImpl(MethodImplOptions.AggressiveOptimization)] public override TSeries Update(TSeries source) { if (source.Count == 0) { return []; } var resultValues = new double[source.Count]; Batch(source.Values, resultValues, LongPeriod, ShortPeriod); var result = new TSeries(); var times = source.Times; for (int i = 0; i < source.Count; i++) { result.Add(new TValue(times[i], resultValues[i])); } // Sync internal state int len = source.Count; if (len >= 2) { var replay = new Pta(LongPeriod, ShortPeriod); for (int i = 0; i < len; i++) { replay.Update(new TValue(times[i], source.Values[i])); } _state = replay._state; } _p_state = _state; return result; } /// Static batch on TSeries. public static TSeries Batch(TSeries source, int longPeriod = 250, int shortPeriod = 40) { var indicator = new Pta(longPeriod, shortPeriod); return indicator.Update(source); } /// /// Static batch calculation on spans. Zero allocation on the hot path. /// [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] public static void Batch(ReadOnlySpan source, Span output, int longPeriod = 250, int shortPeriod = 40) { if (source.Length != output.Length) { throw new ArgumentException("Source and output spans must be of equal length.", nameof(output)); } if (source.Length == 0) { return; } ArgumentOutOfRangeException.ThrowIfLessThan(longPeriod, 3, nameof(longPeriod)); ArgumentOutOfRangeException.ThrowIfLessThan(shortPeriod, 2, nameof(shortPeriod)); if (longPeriod <= shortPeriod) { throw new ArgumentOutOfRangeException(nameof(longPeriod), $"longPeriod ({longPeriod}) must be greater than shortPeriod ({shortPeriod})."); } // Precompute coefficients ComputeHpCoefficients(longPeriod, out double c1L, out double c2L, out double c3L); ComputeHpCoefficients(shortPeriod, out double c1S, out double c2S, out double c3S); // Bar 0 and 1: output = 0 (not enough history for 2nd-order diff) output[0] = 0.0; if (source.Length < 2) { return; } output[1] = 0.0; double hp1 = 0, hp1_1 = 0; double hp2 = 0, hp2_1 = 0; for (int i = 2; i < source.Length; i++) { double diff = source[i] - 2.0 * source[i - 1] + source[i - 2]; double newHp1 = Math.FusedMultiplyAdd(c1L, diff, Math.FusedMultiplyAdd(c2L, hp1, c3L * hp1_1)); double newHp2 = Math.FusedMultiplyAdd(c1S, diff, Math.FusedMultiplyAdd(c2S, hp2, c3S * hp2_1)); output[i] = newHp1 - newHp2; hp1_1 = hp1; hp1 = newHp1; hp2_1 = hp2; hp2 = newHp2; } } /// Calculate factory returning results and indicator. public static (TSeries Results, Pta Indicator) Calculate(TSeries source, int longPeriod = 250, int shortPeriod = 40) { var indicator = new Pta(longPeriod, shortPeriod); TSeries results = indicator.Update(source); return (results, indicator); } public override void Reset() { _state = State.New(); _p_state = _state; } }