Add documentation links for various volatility indicators and channels

- Updated BBWN, BBWP, CCV, CV, CVI, EWMA, GKV, HLV, HV, Jvolty, JVOLTYN, MASSI, NATR, RSV, RV, RVI, TR, UI, VOV, VR, YZV indicators with documentation links.
- Added documentation links for Aberration, Acceleration Bands, Andrews' Pitchfork, Adaptive Price Zone, ATR Bands, Bollinger Bands, Center of Gravity, Donchian Channels, Decay Min-Max Channel, Detrended Synthetic Price, EACP, EBSW, HOMOD, Jurik Volatility Bands, Keltner Channel, MA Envelope, Min-Max Channel, Price Channel, Regression Channels, Standard Deviation Channel, Stoller Average Range Channel, Super Trend Bands, Ultimate Bands, Ultimate Channel, VWAP Bands, and VWAP with Standard Deviation Bands.
This commit is contained in:
Miha Kralj
2026-02-18 11:55:48 -08:00
parent 79c0d72d0a
commit 24e86d762a
332 changed files with 19813 additions and 323 deletions
+132
View File
@@ -0,0 +1,132 @@
using TradingPlatform.BusinessLayer;
namespace QuanTAlib.Tests;
public class RlsIndicatorTests
{
[Fact]
public void RlsIndicator_Constructor_SetsDefaults()
{
var indicator = new RlsIndicator();
Assert.Equal(16, indicator.Order);
Assert.Equal(0.99, indicator.Lambda);
Assert.Equal(SourceType.Close, indicator.Source);
Assert.True(indicator.ShowColdValues);
Assert.Equal("RLS - Recursive Least Squares Adaptive Filter", indicator.Name);
Assert.False(indicator.SeparateWindow);
Assert.True(indicator.OnBackGround);
}
[Fact]
public void RlsIndicator_MinHistoryDepths_EqualsZero()
{
var indicator = new RlsIndicator { Order = 16, Lambda = 0.99 };
Assert.Equal(0, RlsIndicator.MinHistoryDepths);
Assert.Equal(0, ((IWatchlistIndicator)indicator).MinHistoryDepths);
}
[Fact]
public void RlsIndicator_ShortName_IncludesParameters()
{
var indicator = new RlsIndicator { Order = 16, Lambda = 0.99 };
Assert.Contains("RLS", indicator.ShortName, StringComparison.Ordinal);
Assert.Contains("16", indicator.ShortName, StringComparison.Ordinal);
Assert.Contains("0.990", indicator.ShortName, StringComparison.Ordinal);
}
[Fact]
public void RlsIndicator_Initialize_CreatesInternalRls()
{
var indicator = new RlsIndicator { Order = 16, Lambda = 0.99 };
indicator.Initialize();
_ = Assert.Single(indicator.LinesSeries);
}
[Fact]
public void RlsIndicator_ProcessUpdate_HistoricalBar_ComputesValue()
{
var indicator = new RlsIndicator { Order = 4, Lambda = 0.99 };
indicator.Initialize();
var now = DateTime.UtcNow;
indicator.HistoricalData.AddBar(now, 100, 105, 95, 102);
var args = new UpdateArgs(UpdateReason.HistoricalBar);
indicator.ProcessUpdate(args);
Assert.Equal(1, indicator.LinesSeries[0].Count);
Assert.True(double.IsFinite(indicator.LinesSeries[0].GetValue(0)));
}
[Fact]
public void RlsIndicator_ProcessUpdate_NewBar_ComputesValue()
{
var indicator = new RlsIndicator { Order = 4, Lambda = 0.99 };
indicator.Initialize();
var now = DateTime.UtcNow;
indicator.HistoricalData.AddBar(now, 100, 105, 95, 102);
indicator.HistoricalData.AddBar(now.AddMinutes(1), 102, 108, 100, 106);
indicator.ProcessUpdate(new UpdateArgs(UpdateReason.HistoricalBar));
indicator.ProcessUpdate(new UpdateArgs(UpdateReason.NewBar));
Assert.Equal(2, indicator.LinesSeries[0].Count);
}
[Fact]
public void RlsIndicator_ProcessUpdate_NewTick_ProcessesWithoutError()
{
var indicator = new RlsIndicator { Order = 4, Lambda = 0.99 };
indicator.Initialize();
var now = DateTime.UtcNow;
indicator.HistoricalData.AddBar(now, 100, 105, 95, 102);
indicator.ProcessUpdate(new UpdateArgs(UpdateReason.HistoricalBar));
double firstValue = indicator.LinesSeries[0].GetValue(0);
indicator.ProcessUpdate(new UpdateArgs(UpdateReason.NewTick));
double secondValue = indicator.LinesSeries[0].GetValue(0);
Assert.True(double.IsFinite(firstValue));
Assert.True(double.IsFinite(secondValue));
}
[Fact]
public void RlsIndicator_DifferentSourceTypes_Work()
{
var sources = new[] { SourceType.Open, SourceType.High, SourceType.Low, SourceType.Close, SourceType.HL2, SourceType.HLC3 };
foreach (var source in sources)
{
var indicator = new RlsIndicator { Order = 4, Lambda = 0.99, Source = source };
indicator.Initialize();
var now = DateTime.UtcNow;
indicator.HistoricalData.AddBar(now, 100, 110, 90, 105);
indicator.ProcessUpdate(new UpdateArgs(UpdateReason.HistoricalBar));
Assert.True(double.IsFinite(indicator.LinesSeries[0].GetValue(0)),
$"Source {source} should produce finite value");
}
}
[Fact]
public void RlsIndicator_Parameters_CanBeChanged()
{
var indicator = new RlsIndicator { Order = 16, Lambda = 0.99 };
Assert.Equal(16, indicator.Order);
Assert.Equal(0.99, indicator.Lambda);
indicator.Order = 8;
indicator.Lambda = 0.95;
Assert.Equal(8, indicator.Order);
Assert.Equal(0.95, indicator.Lambda);
}
}
+58
View File
@@ -0,0 +1,58 @@
using System.Drawing;
using System.Runtime.CompilerServices;
using TradingPlatform.BusinessLayer;
namespace QuanTAlib;
[SkipLocalsInit]
public sealed class RlsIndicator : Indicator, IWatchlistIndicator
{
[InputParameter("Filter Order (taps)", sortIndex: 1, 2, 64, 1, 0)]
public int Order { get; set; } = 16;
[InputParameter("Forgetting Factor (λ)", sortIndex: 2, 0.9, 1.0, 0.005, 3)]
public double Lambda { get; set; } = 0.99;
[IndicatorExtensions.DataSourceInput]
public SourceType Source { get; set; } = SourceType.Close;
[InputParameter("Show cold values", sortIndex: 21)]
public bool ShowColdValues { get; set; } = true;
private Rls _rls = null!;
private readonly LineSeries _rlsSeries;
private string _sourceName = null!;
private Func<IHistoryItem, double> _priceSelector = null!;
public static int MinHistoryDepths => 0;
int IWatchlistIndicator.MinHistoryDepths => MinHistoryDepths;
public override string ShortName => $"RLS {Order}:{Lambda:F3}:{_sourceName}";
public RlsIndicator()
{
OnBackGround = true;
SeparateWindow = false;
Name = "RLS - Recursive Least Squares Adaptive Filter";
Description = "Adaptive FIR filter with inverse correlation matrix for fast convergence";
_rlsSeries = new LineSeries(name: $"RLS {Order}", color: Color.Orange, width: 2, style: LineStyle.Solid);
AddLineSeries(_rlsSeries);
}
protected override void OnInit()
{
_priceSelector = Source.GetPriceSelector();
_sourceName = Source.ToString();
_rls = new Rls(Order, Lambda);
base.OnInit();
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
protected override void OnUpdate(UpdateArgs args)
{
bool isNew = args.IsNewBar();
var item = HistoricalData[Count - 1, SeekOriginHistory.Begin];
double value = _rls.Update(new TValue(item.TimeLeft.Ticks, _priceSelector(item)), isNew).Value;
_rlsSeries.SetValue(value, _rls.IsHot, ShowColdValues);
}
}
+466
View File
@@ -0,0 +1,466 @@
namespace QuanTAlib;
public class RlsTests
{
private readonly GBM _gbm;
public RlsTests()
{
_gbm = new GBM(startPrice: 100, mu: 0.05, sigma: 0.2, seed: 123);
}
// --- A) Constructor Validation ---
[Fact]
public void Constructor_ValidatesOrder_TooSmall()
{
Assert.Throws<ArgumentOutOfRangeException>(() => new Rls(order: 1));
Assert.Throws<ArgumentOutOfRangeException>(() => new Rls(order: 0));
Assert.Throws<ArgumentOutOfRangeException>(() => new Rls(order: -1));
}
[Fact]
public void Constructor_ValidatesLambda_TooSmall()
{
Assert.Throws<ArgumentOutOfRangeException>(() => new Rls(order: 4, lambda: 0.0));
Assert.Throws<ArgumentOutOfRangeException>(() => new Rls(order: 4, lambda: -0.1));
}
[Fact]
public void Constructor_ValidatesLambda_TooLarge()
{
Assert.Throws<ArgumentOutOfRangeException>(() => new Rls(order: 4, lambda: 1.01));
Assert.Throws<ArgumentOutOfRangeException>(() => new Rls(order: 4, lambda: 2.0));
}
[Fact]
public void Constructor_AcceptsLambdaOne()
{
var ind = new Rls(order: 4, lambda: 1.0);
Assert.Equal(1.0, ind.Lambda);
}
[Fact]
public void Constructor_SetsName()
{
var ind = new Rls(16, 0.99);
Assert.Equal("RLS(16,0.99)", ind.Name);
}
[Fact]
public void Constructor_SetsWarmupPeriod()
{
var ind = new Rls(16, 0.99);
Assert.Equal(17, ind.WarmupPeriod); // order + 1
}
[Fact]
public void Constructor_DefaultParameters()
{
var ind = new Rls();
Assert.Equal(16, ind.Order);
Assert.Equal(0.99, ind.Lambda);
}
[Fact]
public void Constructor_ExposesProperties()
{
var ind = new Rls(8, 0.95);
Assert.Equal(8, ind.Order);
Assert.Equal(0.95, ind.Lambda, 1e-15);
}
// --- B) Basic Calculation ---
[Fact]
public void Calc_ReturnsValue()
{
var ind = new Rls(4, 0.99);
var result = ind.Update(new TValue(DateTime.UtcNow, 100));
Assert.True(double.IsFinite(result.Value));
}
[Fact]
public void Calc_PropertiesAccessible()
{
var ind = new Rls(4, 0.99);
for (int i = 0; i < 10; i++)
{
ind.Update(new TValue(DateTime.UtcNow, 100 + i));
}
Assert.True(double.IsFinite(ind.Last.Value));
Assert.True(ind.IsHot);
Assert.Equal("RLS(4,0.99)", ind.Name);
_ = ind.IsNew;
}
[Fact]
public void Calc_PassthroughDuringWarmup()
{
// During warmup (count <= order), output should equal input
var ind = new Rls(4, 0.99);
for (int i = 0; i < 4; i++)
{
double val = 100 + i;
var result = ind.Update(new TValue(DateTime.UtcNow, val));
Assert.Equal(val, result.Value, 1e-10);
}
}
[Fact]
public void Calc_AdaptiveFilter_FollowsPrice()
{
// RLS is an overlay (price-following) filter — output should track input
var ind = new Rls(8, 0.99);
var data = _gbm.Fetch(200, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1));
double lastInput = 0;
double lastOutput = 0;
foreach (var item in data.Close)
{
lastOutput = ind.Update(item).Value;
lastInput = item.Value;
}
// After adaptation, output should be in the neighborhood of input
double relError = Math.Abs(lastOutput - lastInput) / Math.Abs(lastInput);
Assert.True(relError < 0.5, $"RLS output should track price, relative error = {relError:P2}");
}
// --- C) State + Bar Correction ---
[Fact]
public void Calc_IsNew_AcceptsParameter()
{
// During warmup (passthrough), isNew=false with different value gives different output
var ind = new Rls(4, 0.99);
ind.Update(new TValue(DateTime.UtcNow, 100), isNew: true);
ind.Update(new TValue(DateTime.UtcNow, 105), isNew: true);
double val1 = ind.Last.Value;
// In passthrough mode (count <= order), output = val, so different val = different output
ind.Update(new TValue(DateTime.UtcNow, 110), isNew: false);
double val2 = ind.Last.Value;
Assert.NotEqual(val1, val2);
}
[Fact]
public void Calc_IsNew_False_RollsBackAndRecomputes()
{
// isNew=false should roll back state and recompute with new value
var ind = new Rls(4, 0.99);
for (int i = 0; i < 6; i++)
{
ind.Update(new TValue(DateTime.UtcNow, 100 + i), isNew: true);
}
// Correction with isNew=false
var corrected = ind.Update(new TValue(DateTime.UtcNow, 200), isNew: false);
Assert.True(double.IsFinite(corrected.Value), "Correction should produce finite output");
}
[Fact]
public void IterativeCorrections_RestoreToOriginalState()
{
var ind = new Rls(4, 0.99);
var data = _gbm.Fetch(50, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1));
var series = data.Close;
for (int i = 0; i < series.Count; i++)
{
ind.Update(series[i]);
}
double originalValue = ind.Last.Value;
// Two sequential isNew=false corrections should produce consistent results
var correction1 = ind.Update(new TValue(DateTime.UtcNow, 200), isNew: false);
Assert.True(double.IsFinite(correction1.Value));
var correction2 = ind.Update(new TValue(DateTime.UtcNow, 300), isNew: false);
Assert.True(double.IsFinite(correction2.Value));
// Replaying the same correction value should produce the same result (deterministic)
var correction2b = ind.Update(new TValue(DateTime.UtcNow, 300), isNew: false);
Assert.Equal(correction2.Value, correction2b.Value, 10);
// Replaying original value should restore original prediction
ind.Update(series[^1], isNew: false);
double restoredValue = ind.Last.Value;
Assert.Equal(originalValue, restoredValue, 10);
}
[Fact]
public void Reset_ClearsState()
{
var ind = new Rls(4, 0.99);
var data = _gbm.Fetch(50, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1));
foreach (var item in data.Close)
{
ind.Update(item);
}
ind.Reset();
var ind2 = new Rls(4, 0.99);
var result1 = ind.Update(new TValue(DateTime.UtcNow, 100));
var result2 = ind2.Update(new TValue(DateTime.UtcNow, 100));
Assert.Equal(result2.Value, result1.Value, 10);
}
// --- D) Warmup/Convergence ---
[Fact]
public void IsHot_AfterEnoughBars()
{
var ind = new Rls(4, 0.99);
// Need count > order = 4, so 5 bars
for (int i = 0; i < 4; i++)
{
ind.Update(new TValue(DateTime.UtcNow, 100 + i));
Assert.False(ind.IsHot, $"Should not be hot at count={i + 1}");
}
ind.Update(new TValue(DateTime.UtcNow, 104));
Assert.True(ind.IsHot, "Should be hot after order+1 bars");
// Stays true after more data
for (int i = 0; i < 50; i++)
{
ind.Update(new TValue(DateTime.UtcNow, 100 + i));
}
Assert.True(ind.IsHot);
}
// --- E) Robustness ---
[Fact]
public void NaN_Input_UsesLastValidValue()
{
var ind = new Rls(4, 0.99);
for (int i = 0; i < 6; i++)
{
ind.Update(new TValue(DateTime.UtcNow, 100 + i));
}
var result = ind.Update(new TValue(DateTime.UtcNow, double.NaN));
Assert.True(double.IsFinite(result.Value));
}
[Fact]
public void Infinity_Input_UsesLastValidValue()
{
var ind = new Rls(4, 0.99);
for (int i = 0; i < 6; i++)
{
ind.Update(new TValue(DateTime.UtcNow, 100 + i));
}
var result = ind.Update(new TValue(DateTime.UtcNow, double.PositiveInfinity));
Assert.True(double.IsFinite(result.Value));
var result2 = ind.Update(new TValue(DateTime.UtcNow, double.NegativeInfinity));
Assert.True(double.IsFinite(result2.Value));
}
[Fact]
public void MultipleNaN_ContinuesWithLastValid()
{
var ind = new Rls(4, 0.99);
for (int i = 0; i < 6; i++)
{
ind.Update(new TValue(DateTime.UtcNow, 100 + i));
}
for (int i = 0; i < 10; i++)
{
var result = ind.Update(new TValue(DateTime.UtcNow, double.NaN));
Assert.True(double.IsFinite(result.Value));
}
}
[Fact]
public void BatchCalc_HandlesNaN()
{
double[] input = [100, 105, double.NaN, 110, double.NaN, 115, 120, 125, 130, 135];
double[] output = new double[input.Length];
Rls.Batch(input, output, 4, 0.99);
for (int i = 0; i < output.Length; i++)
{
Assert.True(double.IsFinite(output[i]), $"Output[{i}] should be finite");
}
}
// --- F) Consistency ---
[Fact]
public void AllModes_ProduceSameResult()
{
const int order = 8;
const double lambda = 0.99;
var data = _gbm.Fetch(200, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1));
var series = data.Close;
// 1. Span Mode
double[] spanOutput = new double[series.Count];
Rls.Batch(series.Values.ToArray(), spanOutput, order, lambda);
// 2. TSeries Batch Mode
var batchInd = new Rls(order, lambda);
var batchResult = batchInd.Update(series);
// 3. Streaming Mode
var streamInd = new Rls(order, lambda);
var streamResults = new List<double>();
foreach (var item in series)
{
streamResults.Add(streamInd.Update(item).Value);
}
// 4. Eventing Mode
var pubSource = new TSeries();
var eventInd = new Rls(pubSource, order, lambda);
for (int i = 0; i < series.Count; i++)
{
pubSource.Add(series[i]);
}
// Assert all modes match
for (int i = 0; i < series.Count; i++)
{
Assert.Equal(spanOutput[i], batchResult[i].Value, 1e-9);
Assert.Equal(spanOutput[i], streamResults[i], 1e-9);
}
Assert.Equal(spanOutput[^1], eventInd.Last.Value, 1e-9);
}
// --- G) Span API ---
[Fact]
public void SpanCalc_ValidatesLength()
{
double[] source = new double[10];
double[] output = new double[5]; // Mismatched!
Assert.Throws<ArgumentException>(() => Rls.Batch(source, output));
}
[Fact]
public void SpanCalc_ValidatesOrder()
{
double[] source = new double[10];
double[] output = new double[10];
Assert.Throws<ArgumentOutOfRangeException>(() => Rls.Batch(source, output, order: 1));
}
[Fact]
public void SpanCalc_ValidatesLambda()
{
double[] source = new double[10];
double[] output = new double[10];
Assert.Throws<ArgumentOutOfRangeException>(() => Rls.Batch(source, output, order: 4, lambda: 0.0));
Assert.Throws<ArgumentOutOfRangeException>(() => Rls.Batch(source, output, order: 4, lambda: 1.01));
}
[Fact]
public void SpanCalc_MatchesTSeriesCalc()
{
var data = _gbm.Fetch(200, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1));
var series = data.Close;
// Span
double[] spanOutput = new double[series.Count];
Rls.Batch(series.Values.ToArray(), spanOutput, 8, 0.99);
// TSeries
var ind = new Rls(8, 0.99);
var tseriesResult = ind.Update(series);
for (int i = 0; i < series.Count; i++)
{
Assert.Equal(spanOutput[i], tseriesResult[i].Value, 1e-9);
}
}
[Fact]
public void SpanCalc_NaN_Safe()
{
double[] input = new double[50];
for (int i = 0; i < 50; i++)
{
input[i] = i % 7 == 0 ? double.NaN : 100.0 + Math.Sin(i * 0.1);
}
double[] output = new double[50];
Rls.Batch(input, output, 4, 0.99);
for (int i = 0; i < output.Length; i++)
{
Assert.True(double.IsFinite(output[i]), $"Output[{i}] should be finite with NaN input");
}
}
// --- H) Chainability ---
[Fact]
public void Pub_FiresOnUpdate()
{
var ind = new Rls(4, 0.99);
int fireCount = 0;
ind.Pub += (object? _, in TValueEventArgs _) => fireCount++;
ind.Update(new TValue(DateTime.UtcNow, 100));
ind.Update(new TValue(DateTime.UtcNow, 105));
Assert.Equal(2, fireCount);
}
[Fact]
public void EventChaining_Works()
{
var source = new TSeries();
var ind = new Rls(source, 4, 0.99);
source.Add(new TValue(DateTime.UtcNow, 100));
source.Add(new TValue(DateTime.UtcNow, 105));
Assert.True(double.IsFinite(ind.Last.Value));
}
// --- Additional ---
[Fact]
public void DifferentParameters_ProduceDifferentResults()
{
var data = _gbm.Fetch(200, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1));
var series = data.Close;
var ind1 = new Rls(8, 0.99);
var ind2 = new Rls(16, 0.95);
foreach (var item in series)
{
ind1.Update(item);
ind2.Update(item);
}
Assert.NotEqual(ind1.Last.Value, ind2.Last.Value);
}
[Fact]
public void LargeDataset_DoesNotThrow()
{
var data = _gbm.Fetch(5000, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1));
var series = data.Close;
double[] input = series.Values.ToArray();
double[] output = new double[input.Length];
Rls.Batch(input, output, 16, 0.99);
Assert.True(double.IsFinite(output[^1]));
}
}
+228
View File
@@ -0,0 +1,228 @@
namespace QuanTAlib.Tests;
/// <summary>
/// Validation tests for the RLS Adaptive Filter.
/// Since RLS is a custom adaptive filter with no direct external library equivalent,
/// validation uses self-consistency: adaptive convergence, streaming/span parity,
/// determinism, stability, and mathematical properties of the RLS algorithm.
/// RLS should converge faster than LMS due to the inverse correlation matrix.
/// </summary>
public class RlsValidationTests
{
[Fact]
public void Validate_AdaptiveConvergence_SineWave()
{
// RLS should learn to predict a periodic signal with decreasing error
const int T = 500;
double[] sine = new double[T];
for (int i = 0; i < T; i++)
{
sine[i] = 100.0 + 10.0 * Math.Sin(2 * Math.PI * i / 40.0);
}
double[] output = new double[T];
Rls.Batch(sine, output, 8, 0.99);
// Compute mean squared error in first quarter vs last quarter
double mseFirst = 0, mseLast = 0;
int q = T / 4;
for (int i = 0; i < q; i++)
{
double e = sine[i] - output[i];
mseFirst += e * e;
}
for (int i = T - q; i < T; i++)
{
double e = sine[i] - output[i];
mseLast += e * e;
}
mseFirst /= q;
mseLast /= q;
Assert.True(mseLast < mseFirst, $"Error should decrease: first quarter MSE={mseFirst:F4}, last quarter MSE={mseLast:F4}");
}
[Fact]
public void Validate_StreamingMatchesSpan()
{
var gbm = new GBM(startPrice: 100, mu: 0.05, sigma: 0.2, seed: 42);
var data = gbm.Fetch(500, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1));
double[] input = data.Close.Values.ToArray();
// Span path
double[] spanOut = new double[input.Length];
Rls.Batch(input, spanOut, 8, 0.99);
// Streaming path
var ind = new Rls(8, 0.99);
var streamResults = new double[input.Length];
for (int i = 0; i < input.Length; i++)
{
streamResults[i] = ind.Update(new TValue(DateTime.UtcNow, input[i])).Value;
}
for (int i = 0; i < input.Length; i++)
{
Assert.Equal(spanOut[i], streamResults[i], 1e-9);
}
}
[Fact]
public void Validate_ConstantInput_ConvergesToConstant()
{
// Constant input -> filter should predict constant -> output ~ input after warmup
double[] input = Enumerable.Repeat(50.0, 500).ToArray();
double[] output = new double[500];
Rls.Batch(input, output, 8, 0.99);
// After warmup, output should converge close to input
Assert.True(Math.Abs(output[^1] - 50.0) < 1.0,
$"Constant input should yield ~50, got {output[^1]}");
}
[Fact]
public void Validate_Deterministic()
{
var gbm = new GBM(startPrice: 100, mu: 0.05, sigma: 0.2, seed: 99);
var data = gbm.Fetch(500, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1));
double[] input = data.Close.Values.ToArray();
double[] out1 = new double[input.Length];
double[] out2 = new double[input.Length];
Rls.Batch(input, out1, 8, 0.99);
Rls.Batch(input, out2, 8, 0.99);
for (int i = 0; i < input.Length; i++)
{
Assert.Equal(out1[i], out2[i], 15);
}
}
[Fact]
public void Validate_OutputFollowsInput()
{
// RLS is an overlay filter — output should track input direction
var gbm = new GBM(startPrice: 100, mu: 0.05, sigma: 0.2, seed: 77);
var data = gbm.Fetch(1000, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1));
double[] input = data.Close.Values.ToArray();
double[] output = new double[input.Length];
Rls.Batch(input, output, 8, 0.99);
// Correlation between input and output should be positive and high
double meanIn = 0, meanOut = 0;
int start = 50; // skip warmup
int n = input.Length - start;
for (int i = start; i < input.Length; i++)
{
meanIn += input[i];
meanOut += output[i];
}
meanIn /= n;
meanOut /= n;
double cov = 0, varIn = 0, varOut = 0;
for (int i = start; i < input.Length; i++)
{
double dIn = input[i] - meanIn;
double dOut = output[i] - meanOut;
cov += dIn * dOut;
varIn += dIn * dIn;
varOut += dOut * dOut;
}
double corr = cov / Math.Sqrt(varIn * varOut);
Assert.True(corr > 0.5, $"Output should track input, correlation = {corr:F4}");
}
[Fact]
public void Validate_LargeDataset_Stable()
{
var gbm = new GBM(startPrice: 100, mu: 0.05, sigma: 0.2, seed: 55);
var data = gbm.Fetch(5000, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1));
double[] input = data.Close.Values.ToArray();
double[] output = new double[input.Length];
Rls.Batch(input, output, 8, 0.99);
for (int i = 0; i < output.Length; i++)
{
Assert.True(double.IsFinite(output[i]), $"Output[{i}] is not finite: {output[i]}");
}
}
[Fact]
public void Validate_NaN_Batch_Safe()
{
double[] input = new double[100];
for (int i = 0; i < 100; i++)
{
input[i] = i % 7 == 0 ? double.NaN : 100.0 + Math.Sin(i * 0.1);
}
double[] output = new double[100];
Rls.Batch(input, output, 4, 0.99);
for (int i = 0; i < output.Length; i++)
{
Assert.True(double.IsFinite(output[i]), $"Output[{i}] should be finite with NaN input");
}
}
[Fact]
public void Validate_FasterConvergenceThanLMS()
{
// RLS should converge faster than LMS on a step function
const int T = 200;
double[] input = new double[T];
for (int i = 0; i < T; i++)
{
input[i] = i < 50 ? 100.0 : 120.0;
}
double[] rlsOut = new double[T];
double[] lmsOut = new double[T];
Rls.Batch(input, rlsOut, 4, 0.99);
Lms.Batch(input, lmsOut, 4, 0.5);
// Measure error in the adaptation window (bars 55-70 after step)
double rlsErr = 0, lmsErr = 0;
for (int i = 55; i < 70; i++)
{
rlsErr += Math.Abs(rlsOut[i] - 120.0);
lmsErr += Math.Abs(lmsOut[i] - 120.0);
}
Assert.True(rlsErr < lmsErr,
$"RLS should converge faster: RLS error={rlsErr:F4}, LMS error={lmsErr:F4}");
}
[Fact]
public void Validate_DifferentLambda_ProduceDifferentOutput()
{
var gbm = new GBM(startPrice: 100, mu: 0.05, sigma: 0.2, seed: 33);
var data = gbm.Fetch(500, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1));
double[] input = data.Close.Values.ToArray();
double[] out1 = new double[input.Length];
double[] out2 = new double[input.Length];
Rls.Batch(input, out1, 8, 0.99);
Rls.Batch(input, out2, 8, 0.95);
bool anyDifferent = false;
for (int i = 20; i < input.Length; i++)
{
if (Math.Abs(out1[i] - out2[i]) > 1e-12)
{
anyDifferent = true;
break;
}
}
Assert.True(anyDifferent, "Different lambda values should produce different output");
}
}
+449
View File
@@ -0,0 +1,449 @@
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
namespace QuanTAlib;
/// <summary>
/// RLS: Recursive Least Squares Adaptive Filter
/// An adaptive FIR filter that maintains an inverse correlation matrix P to achieve
/// faster convergence than LMS. Uses a forgetting factor λ to control memory horizon.
/// Converges in ~order iterations with O(order²) per-bar complexity.
/// </summary>
/// <remarks>
/// The algorithm is based on a Pine Script implementation:
/// https://github.com/mihakralj/pinescript/blob/main/indicators/filters/rls.md
///
/// Key properties:
/// - Adaptive FIR: weight vector w[0..order-1] learns from streaming data
/// - Predicts src[0] from src[1]..src[order] (no look-ahead)
/// - Gain vector: k = P·x / (λ + x^T·P·x)
/// - P update: P = (1/λ)(P - k·(P·x)^T)
/// - Overlay indicator (price-following)
/// - O(order²) per bar for both prediction and weight/matrix update
///
/// Complexity: O(order²) per bar
/// </remarks>
[SkipLocalsInit]
public sealed class Rls : AbstractBase
{
private const double Epsilon = 1e-30;
private readonly int _order;
private readonly double _lambda;
private readonly double _invLambda;
private readonly RingBuffer _inputBuffer;
private readonly double[] _weights;
private readonly double[] _p_weights;
private readonly double[] _P; // inverse correlation matrix (order x order), row-major
private readonly double[] _p_P; // snapshot for bar correction
private ITValuePublisher? _publisher;
private TValuePublishedHandler? _handler;
private bool _isNew;
[StructLayout(LayoutKind.Auto)]
private record struct State
{
public double LastValid;
public int Count;
}
private State _state;
private State _p_state;
/// <summary>Number of FIR taps (adaptive weights).</summary>
public int Order => _order;
/// <summary>Forgetting factor controlling memory horizon (0 &lt; λ ≤ 1).</summary>
public double Lambda => _lambda;
public bool IsNew => _isNew;
public override bool IsHot => _state.Count > _order;
public Rls(int order = 16, double lambda = 0.99)
{
if (order < 2)
{
throw new ArgumentOutOfRangeException(nameof(order), "Filter order must be >= 2.");
}
if (lambda <= 0.0 || lambda > 1.0)
{
throw new ArgumentOutOfRangeException(nameof(lambda), "Forgetting factor lambda must be in (0, 1].");
}
_order = order;
_lambda = lambda;
_invLambda = 1.0 / lambda;
Name = $"RLS({order},{lambda:F2})";
WarmupPeriod = order + 1;
// Weight vector + snapshot for bar correction
_weights = new double[order];
_p_weights = new double[order];
// Inverse correlation matrix P = delta * I (high initial uncertainty)
const double delta = 100.0;
int matSize = order * order;
_P = new double[matSize];
_p_P = new double[matSize];
for (int i = 0; i < order; i++)
{
_P[i * order + i] = delta;
}
// Ring buffer holds order+1 values: current + order past values
_inputBuffer = new RingBuffer(order + 1);
_state.LastValid = double.NaN;
}
public Rls(ITValuePublisher source, int order = 16, double lambda = 0.99)
: this(order, lambda)
{
_publisher = source;
_handler = Handle;
source.Pub += _handler;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private void Handle(object? sender, in TValueEventArgs args)
{
Update(args.Value, args.IsNew);
}
public override TSeries Update(TSeries source)
{
if (source.Count == 0)
{
return [];
}
double[] values = source.Values.ToArray();
double[] results = new double[values.Length];
Batch(values, results, _order, _lambda);
TSeries output = [];
for (int i = 0; i < values.Length; i++)
{
output.Add(source[i].Time, results[i]);
}
// Resync internal state by replaying
Reset();
for (int i = 0; i < source.Count; i++)
{
Update(source[i]);
}
return output;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public override TValue Update(TValue input, bool isNew = true)
{
_isNew = isNew;
int matSize = _order * _order;
if (isNew)
{
_p_state = _state;
Array.Copy(_weights, _p_weights, _order);
Array.Copy(_P, _p_P, matSize);
}
else
{
_state = _p_state;
Array.Copy(_p_weights, _weights, _order);
Array.Copy(_p_P, _P, matSize);
}
var s = _state;
// Handle bad data — last-valid substitution
double val = input.Value;
if (!double.IsFinite(val))
{
val = double.IsFinite(s.LastValid) ? s.LastValid : 0.0;
}
else
{
s.LastValid = val;
}
// Input buffer: Add for new bars, UpdateNewest for corrections
if (isNew)
{
_inputBuffer.Add(val);
}
else
{
_inputBuffer.UpdateNewest(val);
}
double result;
if (_inputBuffer.Count <= _order)
{
// Not enough history to form prediction — pass through
result = val;
}
else
{
// --- Step 1: Prediction y = w^T * x ---
double y = 0.0;
for (int i = 0; i < _order; i++)
{
double xi = _inputBuffer[^(i + 2)]; // src[i+1]
y = Math.FusedMultiplyAdd(_weights[i], xi, y);
}
// --- RLS update — only learn from confirmed bars ---
if (isNew)
{
// --- Step 2: Compute Px = P * x ---
// skipcq: CS-W1082 - stackalloc safe: order is bounded by constructor validation
Span<double> px = stackalloc double[_order];
for (int i = 0; i < _order; i++)
{
double rowSum = 0.0;
int rowBase = i * _order;
for (int j = 0; j < _order; j++)
{
double xj = _inputBuffer[^(j + 2)];
rowSum = Math.FusedMultiplyAdd(_P[rowBase + j], xj, rowSum);
}
px[i] = rowSum;
}
// --- Step 3: Compute denom = λ + x^T * Px ---
double denom = _lambda;
for (int i = 0; i < _order; i++)
{
double xi = _inputBuffer[^(i + 2)];
denom = Math.FusedMultiplyAdd(xi, px[i], denom);
}
// --- Step 4: Gain vector k = Px / denom ---
double invDenom = denom > Epsilon ? 1.0 / denom : 0.0;
// skipcq: CS-W1082 - stackalloc safe: order is bounded
Span<double> k = stackalloc double[_order];
for (int i = 0; i < _order; i++)
{
k[i] = px[i] * invDenom;
}
// --- Step 5: Weight update w = w + k * error ---
double error = val - y;
for (int i = 0; i < _order; i++)
{
_weights[i] = Math.FusedMultiplyAdd(k[i], error, _weights[i]);
}
// --- Step 6: P update: P = (1/λ)(P - k * Px^T) ---
for (int i = 0; i < _order; i++)
{
double ki = k[i];
int rowBase = i * _order;
for (int j = 0; j < _order; j++)
{
_P[rowBase + j] = _invLambda * (_P[rowBase + j] - ki * px[j]);
}
}
}
result = y;
}
if (isNew)
{
s.Count++;
}
_state = s;
Last = new TValue(input.Time, result);
PubEvent(Last, isNew);
return Last;
}
public static TSeries Batch(TSeries source, int order = 16, double lambda = 0.99)
{
var indicator = new Rls(order, lambda);
return indicator.Update(source);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void Batch(ReadOnlySpan<double> source, Span<double> output,
int order = 16, double lambda = 0.99)
{
if (source.Length != output.Length)
{
throw new ArgumentException("Source and output spans must be of the same length.", nameof(output));
}
if (order < 2)
{
throw new ArgumentOutOfRangeException(nameof(order), "Filter order must be >= 2.");
}
if (lambda <= 0.0 || lambda > 1.0)
{
throw new ArgumentOutOfRangeException(nameof(lambda), "Forgetting factor lambda must be in (0, 1].");
}
double invLambda = 1.0 / lambda;
// Weight vector
double[] w = new double[order];
var ring = new RingBuffer(order + 1);
// Inverse correlation matrix P = delta * I
const double delta = 100.0;
int matSize = order * order;
double[] P = new double[matSize];
for (int i = 0; i < order; i++)
{
P[i * order + i] = delta;
}
// Temporary buffers for Px and k
double[] px = new double[order];
double[] k = new double[order];
double lastValid = 0;
if (source.Length > 0)
{
lastValid = source[0];
if (!double.IsFinite(lastValid))
{
lastValid = 0;
}
}
for (int n = 0; n < source.Length; n++)
{
double val = source[n];
if (!double.IsFinite(val))
{
val = lastValid;
}
else
{
lastValid = val;
}
ring.Add(val, true);
if (ring.Count <= order)
{
output[n] = val;
continue;
}
// Step 1: Prediction y = w^T * x
double y = 0.0;
for (int i = 0; i < order; i++)
{
double xi = ring[^(i + 2)];
y = Math.FusedMultiplyAdd(w[i], xi, y);
}
// Step 2: Px = P * x
for (int i = 0; i < order; i++)
{
double rowSum = 0.0;
int rowBase = i * order;
for (int j = 0; j < order; j++)
{
double xj = ring[^(j + 2)];
rowSum = Math.FusedMultiplyAdd(P[rowBase + j], xj, rowSum);
}
px[i] = rowSum;
}
// Step 3: denom = λ + x^T * Px
double denom = lambda;
for (int i = 0; i < order; i++)
{
double xi = ring[^(i + 2)];
denom = Math.FusedMultiplyAdd(xi, px[i], denom);
}
// Step 4: k = Px / denom
double invDenom = denom > Epsilon ? 1.0 / denom : 0.0;
for (int i = 0; i < order; i++)
{
k[i] = px[i] * invDenom;
}
// Step 5: w = w + k * error
double error = val - y;
for (int i = 0; i < order; i++)
{
w[i] = Math.FusedMultiplyAdd(k[i], error, w[i]);
}
// Step 6: P = (1/λ)(P - k * Px^T)
for (int i = 0; i < order; i++)
{
double ki = k[i];
int rowBase = i * order;
for (int j = 0; j < order; j++)
{
P[rowBase + j] = invLambda * (P[rowBase + j] - ki * px[j]);
}
}
output[n] = y;
}
}
public override void Reset()
{
_state = default;
_state.LastValid = double.NaN;
_p_state = default;
_inputBuffer.Clear();
Array.Clear(_weights);
Array.Clear(_p_weights);
// Reset P to delta * I
const double delta = 100.0;
Array.Clear(_P);
Array.Clear(_p_P);
for (int i = 0; i < _order; i++)
{
_P[i * _order + i] = delta;
}
Last = default;
}
public override void Prime(ReadOnlySpan<double> source, TimeSpan? step = null)
{
foreach (double val in source)
{
Update(new TValue(DateTime.UtcNow, val), isNew: true);
}
}
public static (TSeries Results, Rls Indicator) Calculate(TSeries source,
int order = 16, double lambda = 0.99)
{
var indicator = new Rls(order, lambda);
TSeries results = indicator.Update(source);
return (results, indicator);
}
protected override void Dispose(bool disposing)
{
if (disposing && _publisher != null && _handler != null)
{
_publisher.Pub -= _handler;
_publisher = null;
_handler = null;
}
base.Dispose(disposing);
}
}
+154
View File
@@ -0,0 +1,154 @@
# RLS: Recursive Least Squares Adaptive Filter
> "The man who has no patience has no wisdom." — but waiting is not the same as convergence. RLS converges where LMS merely approaches.
## Introduction
The Recursive Least Squares (RLS) adaptive filter is the Rolls-Royce of adaptive FIR filters. Where LMS crawls toward the Wiener solution one gradient step at a time, RLS arrives in approximately *order* iterations by maintaining an inverse correlation matrix $P$ that captures the full second-order statistics of the input signal. The trade-off is computational: $O(n^2)$ per bar versus LMS's $O(n)$, where $n$ is the filter order. For orders below 64, the convergence advantage typically outweighs the cost.
## Historical Context
RLS traces its lineage to Gauss's method of least squares (1795) and Kalman's recursive state estimation (1960). The exponentially-weighted RLS form — with forgetting factor $\lambda$ — was formalized in the signal processing literature of the 1970s and 1980s, primarily by Haykin, Widrow, and Ljung. Unlike LMS, which adapts proportionally to the instantaneous gradient, RLS minimizes the weighted sum of all past squared errors, making it optimal in a least-squares sense at every time step.
In financial applications, RLS excels at tracking non-stationary price dynamics. The forgetting factor $\lambda$ controls the effective memory horizon: $\lambda = 0.99$ gives a memory of roughly $1/(1-\lambda) = 100$ bars, while $\lambda = 0.95$ compresses memory to 20 bars. This makes RLS particularly suited for regime changes and structural breaks where LMS's fixed step size is too slow to react.
The implementation here follows the standard RLS algorithm with no look-ahead, no leakage, and no regularization beyond the initial $P = \delta I$ scaling.
## Architecture and Physics
### 1. Adaptive Weight Vector
The filter maintains $n$ weights $w_0, w_1, \ldots, w_{n-1}$ that adapt to predict the current input from its recent history:
$$\hat{y}(t) = \sum_{i=0}^{n-1} w_i \cdot x(t-i-1)$$
The prediction uses values $x(t-1)$ through $x(t-n)$ — no look-ahead.
### 2. Inverse Correlation Matrix
The core of RLS is the $n \times n$ inverse correlation matrix $P$, initialized to $\delta I$ where $\delta = 100$ represents high initial uncertainty. This matrix is updated recursively at each step, avoiding the $O(n^3)$ cost of explicit matrix inversion.
### 3. Gain Vector
The Kalman-like gain vector determines how much each weight adjusts in response to prediction error:
$$k(t) = \frac{P(t-1) \cdot x(t)}{\lambda + x(t)^T \cdot P(t-1) \cdot x(t)}$$
### 4. Weight and Matrix Update
After computing the a priori error $e(t) = d(t) - \hat{y}(t)$:
$$w(t) = w(t-1) + k(t) \cdot e(t)$$
$$P(t) = \frac{1}{\lambda}\left(P(t-1) - k(t) \cdot x(t)^T \cdot P(t-1)\right)$$
### 5. Forgetting Factor
The forgetting factor $\lambda \in (0, 1]$ exponentially discounts past observations. The effective memory window is approximately $1/(1-\lambda)$ samples.
| $\lambda$ | Effective Memory | Use Case |
|-----------|------------------|----------|
| 1.00 | Infinite (growing) | Stationary signals |
| 0.99 | ~100 bars | Moderate non-stationarity |
| 0.95 | ~20 bars | Fast-changing dynamics |
| 0.90 | ~10 bars | Highly non-stationary |
## Mathematical Foundation
### Transfer Function (z-domain)
RLS is a time-varying FIR filter. At convergence on a stationary signal, the weight vector approaches the Wiener solution:
$$w_{opt} = R^{-1} p$$
where $R$ is the input autocorrelation matrix and $p$ is the cross-correlation vector between input and desired signal. The z-domain transfer function at convergence is:
$$H(z) = \sum_{i=0}^{n-1} w_i \cdot z^{-(i+1)}$$
### Convergence Analysis
RLS converges in approximately $n$ iterations (where $n$ is the filter order), compared to LMS which requires $O(n / \mu_{\text{eff}})$ iterations. This is because RLS effectively pre-whitens the input through the $P$ matrix, decorrelating the gradient components.
### Stability Condition
The algorithm is stable when $0 < \lambda \leq 1$ and $\delta > 0$. The initial $P = \delta I$ determines convergence speed: larger $\delta$ means faster initial adaptation but potentially larger transient errors.
### Parameter Mapping
| Parameter | Symbol | Default | Range | Effect |
|-----------|--------|---------|-------|--------|
| Order | $n$ | 16 | $[2, 64]$ | Filter taps; higher = more modeling capacity |
| Lambda | $\lambda$ | 0.99 | $(0, 1]$ | Forgetting factor; lower = shorter memory |
| Delta | $\delta$ | 100.0 | $(0, \infty)$ | Initial P scaling; higher = faster initial adaptation |
## Performance Profile
### Operation Count Per Bar
| Operation | Count | Notes |
|-----------|-------|-------|
| Prediction ($w^T x$) | $O(n)$ | FMA inner product |
| $P \cdot x$ | $O(n^2)$ | Matrix-vector multiply |
| Gain vector $k$ | $O(n)$ | Scalar division + scale |
| Weight update | $O(n)$ | $w += k \cdot e$ |
| P update | $O(n^2)$ | Rank-1 outer product subtraction |
| **Total** | **$O(n^2)$** | Dominated by P operations |
### Memory Usage
| Component | Size | Notes |
|-----------|------|-------|
| Weights $w$ | $2n$ doubles | Current + snapshot |
| Matrix $P$ | $2n^2$ doubles | Current + snapshot |
| Input buffer | $n+1$ doubles | RingBuffer |
| **Total** | **$2n^2 + 3n + 1$** | ~4 KB for order=16 |
### Quality Metrics
| Metric | Score (1-10) | Notes |
|--------|:---:|-------|
| Smoothness | 7 | Tracks signal closely |
| Lag | 2 | Minimal prediction lag |
| Overshoot | 4 | Can overshoot in transients |
| Noise rejection | 7 | Good with appropriate $\lambda$ |
| Adaptability | 9 | Fast convergence to optimal |
| Computational cost | 4 | $O(n^2)$ limits practical order |
## Validation
RLS is a custom adaptive filter with no direct equivalent in standard TA libraries. Validation uses self-consistency tests.
| Test | Method | Result |
|------|--------|--------|
| Convergence | MSE decreases on sine wave | First quarter MSE > last quarter MSE |
| Streaming = Span | Mode parity | Match to $10^{-9}$ |
| Determinism | Two identical runs | Match to $10^{-15}$ |
| Constant input | Converge to constant | $\|{y - 50}\| < 1$ |
| Price tracking | Correlation test | $r > 0.5$ |
| Stability | 5000-bar dataset | All outputs finite |
| NaN safety | Interspersed NaN | All outputs finite |
| Faster than LMS | Step response comparison | RLS error < LMS error |
| Lambda sensitivity | Different $\lambda$ values | Different outputs |
## Common Pitfalls
1. **Order too large.** RLS is $O(n^2)$; order 64 means 4096 multiply-adds per bar for the P update alone. Keep order ≤ 32 for real-time use. Impact: 4× latency per doubling of order.
2. **Lambda too small.** Values below 0.9 create a memory horizon of fewer than 10 bars, causing wild weight oscillations. The filter "forgets" useful history and tracks noise. Impact: output variance increases by 3-5×.
3. **P matrix blowup.** Without regularization, $P$ can grow unbounded when $\lambda < 1$ and the input lacks sufficient excitation. The implementation guards against this via the $\epsilon$-denominator clamp. Impact: numerical overflow → NaN propagation.
4. **Confusing lambda with LMS mu.** Lambda is a forgetting factor (higher = more memory), while LMS mu is a step size (higher = faster adaptation). They have opposite semantics despite both controlling adaptation speed.
5. **Initial transient.** The first ~order bars produce passthrough output while the buffer fills. The P matrix starts at $\delta I$, so the first few predictions after warmup may be large. Impact: 2-5 bars of unreliable output after warmup.
6. **Bar correction cost.** Each isNew=false correction requires restoring both the weight vector ($n$ copies) and the P matrix ($n^2$ copies). For order=16, that is 256+16 = 272 doubles copied per correction. Impact: correction cost proportional to $n^2$.
7. **Not suitable for SIMD.** The sequential dependency chain (P update depends on gain, gain depends on P·x) prevents vectorization of the inner loop. Unlike simple FIR filters, RLS cannot benefit from AVX2/SSE parallelism.
## References
- Haykin, S. (2002). *Adaptive Filter Theory*. 4th ed. Prentice Hall. Chapters 9-10.
- Ljung, L. & Soderstrom, T. (1983). *Theory and Practice of Recursive Identification*. MIT Press.
- Sayed, A.H. (2008). *Adaptive Filters*. Wiley-IEEE Press.
- Kalman, R.E. (1960). "A New Approach to Linear Filtering and Prediction Problems." *Journal of Basic Engineering*, 82(1), 35-45.
+96
View File
@@ -0,0 +1,96 @@
// The MIT License (MIT)
// © mihakralj
//@version=6
indicator("Recursive Least Squares Adaptive Filter (RLS)", "RLS", overlay=true)
//@function Applies RLS adaptive FIR filter to input series
//@param src Input series to filter
//@param order Number of FIR filter taps (adaptive weights)
//@param lambda Forgetting factor (0 < lambda <= 1). Controls memory horizon.
//@returns Adaptively filtered series
//@optimized Uses inverse correlation matrix for O(order²) convergence per bar,
// faster convergence than LMS at higher computational cost
rls(series float src, simple int order, simple float lambda) =>
if order < 2
runtime.error("Filter order must be >= 2")
if lambda <= 0.0 or lambda > 1.0
runtime.error("Forgetting factor lambda must be in (0, 1]")
// weight vector w[order]
var array<float> w = array.new_float(order, 0.0)
// inverse correlation matrix P[order x order], stored row-major
// initialized to delta * I (large diagonal = high initial uncertainty)
var float delta = 100.0
var array<float> P = array.new_float(order * order, 0.0)
var bool initialized = false
if not initialized
for i = 0 to order - 1
array.set(P, i * order + i, delta)
initialized := true
// build input vector x = [src[1], src[2], ..., src[order]]
array<float> x = array.new_float(order, 0.0)
for i = 0 to order - 1
array.set(x, i, nz(src[i + 1], 0.0))
// --- Step 1: prediction ---
// y = w^T * x
float y = 0.0
for i = 0 to order - 1
y += array.get(w, i) * array.get(x, i)
// --- Step 2: a priori error ---
float e = nz(src, 0.0) - y
// --- Step 3: gain vector k = P*x / (lambda + x^T*P*x) ---
// compute Px = P * x
array<float> Px = array.new_float(order, 0.0)
for i = 0 to order - 1
float row_sum = 0.0
for j = 0 to order - 1
row_sum += array.get(P, i * order + j) * array.get(x, j)
array.set(Px, i, row_sum)
// compute denom = lambda + x^T * Px
float denom = lambda
for i = 0 to order - 1
denom += array.get(x, i) * array.get(Px, i)
// k = Px / denom
array<float> k = array.new_float(order, 0.0)
float inv_denom = denom > 1e-30 ? 1.0 / denom : 0.0
for i = 0 to order - 1
array.set(k, i, array.get(Px, i) * inv_denom)
// --- Step 4: weight update ---
// w = w + k * e
for i = 0 to order - 1
array.set(w, i, array.get(w, i) + array.get(k, i) * e)
// --- Step 5: P update ---
// P = (1/lambda) * (P - k * x^T * P)
// equivalently: P = (1/lambda) * (P - k * Px^T)
// since Px = P*x, then k*x^T*P = k*(Px)^T (row-by-row outer product)
float inv_lambda = 1.0 / lambda
for i = 0 to order - 1
float ki = array.get(k, i)
for j = 0 to order - 1
float old_pij = array.get(P, i * order + j)
float new_pij = inv_lambda * (old_pij - ki * array.get(Px, j))
array.set(P, i * order + j, new_pij)
y
// ---------- Main loop ----------
// Inputs
i_order = input.int(16, "Filter Order (taps)", minval=2, maxval=64)
i_lambda = input.float(0.99, "Forgetting Factor (λ)", minval=0.9, maxval=1.0, step=0.005)
i_source = input.source(close, "Source")
// Calculation
rls_val = rls(i_source, i_order, i_lambda)
// Plot
plot(rls_val, "RLS", color=color.orange, linewidth=2)