mirror of
https://github.com/mihakralj/QuanTAlib.git
synced 2026-08-05 04:27:43 +00:00
1013 lines
31 KiB
C#
1013 lines
31 KiB
C#
|
using System.Buffers;
|
|||
|
|
using System.Runtime.CompilerServices;
|
|||
|
|
using System.Runtime.InteropServices;
|
|||
|
|
|
|||
|
|
namespace QuanTAlib;
|
|||
|
|
|
|||
|
|
/// <summary>
|
|||
|
|
/// ADF: Augmented Dickey-Fuller unit root test.
|
|||
|
|
/// </summary>
|
|||
|
|
/// <remarks>
|
|||
|
|
/// Tests the null hypothesis that a time series has a unit root (is non-stationary).
|
|||
|
|
/// Output is the MacKinnon (1994) approximate p-value: values near 0 indicate
|
|||
|
|
/// stationarity, values near 1 indicate a unit root (non-stationary / random walk).
|
|||
|
|
///
|
|||
|
|
/// The full regression model is:
|
|||
|
|
/// Δy_t = α + β·t + γ·y_{t-1} + Σ(δ_i·Δy_{t-i}) + ε_t
|
|||
|
|
///
|
|||
|
|
/// where γ is the coefficient of interest. The test statistic is t = γ̂ / SE(γ̂).
|
|||
|
|
/// P-value is computed via MacKinnon (1994) polynomial interpolation with standard
|
|||
|
|
/// normal CDF.
|
|||
|
|
///
|
|||
|
|
/// Complexity: O(period × maxLag²) per update — OLS solve is O(k³) with k ≤ 8.
|
|||
|
|
/// </remarks>
|
|||
|
|
[SkipLocalsInit]
|
|||
|
|
public sealed class Adf : AbstractBase
|
|||
|
|
{
|
|||
|
|
/// <summary>
|
|||
|
|
/// Regression model type for the ADF test.
|
|||
|
|
/// </summary>
|
|||
|
|
public enum AdfRegression
|
|||
|
|
{
|
|||
|
|
/// <summary>No constant, no trend: Δy_t = γ·y_{t-1} + Σ(δ_i·Δy_{t-i}) + ε_t</summary>
|
|||
|
|
NoConstant = 0,
|
|||
|
|
/// <summary>Constant only: Δy_t = α + γ·y_{t-1} + Σ(δ_i·Δy_{t-i}) + ε_t</summary>
|
|||
|
|
Constant = 1,
|
|||
|
|
/// <summary>Constant and trend: Δy_t = α + β·t + γ·y_{t-1} + Σ(δ_i·Δy_{t-i}) + ε_t</summary>
|
|||
|
|
ConstantAndTrend = 2
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
private const double InvSqrt2 = 0.70710678118654752; // 1/√2
|
|||
|
|
private const int MinPeriod = 20;
|
|||
|
|
private const int MaxRegressors = 8; // max k for stackalloc Cholesky
|
|||
|
|
|
|||
|
|
private readonly int _period;
|
|||
|
|
private readonly int _maxLag;
|
|||
|
|
private readonly AdfRegression _regression;
|
|||
|
|
private readonly RingBuffer _buffer;
|
|||
|
|
private double _lastValidValue;
|
|||
|
|
private int _inputCount;
|
|||
|
|
private int _inputCountSaved;
|
|||
|
|
|
|||
|
|
/// <summary>Test statistic (t-value for γ̂).</summary>
|
|||
|
|
public double Statistic { get; private set; }
|
|||
|
|
|
|||
|
|
/// <summary>P-value from MacKinnon approximation.</summary>
|
|||
|
|
public double PValue { get; private set; }
|
|||
|
|
|
|||
|
|
/// <summary>Number of lags used in the augmented regression.</summary>
|
|||
|
|
public int LagsUsed { get; private set; }
|
|||
|
|
|
|||
|
|
public override bool IsHot => _inputCount >= _period;
|
|||
|
|
|
|||
|
|
/// <summary>
|
|||
|
|
/// Creates a new Augmented Dickey-Fuller test indicator.
|
|||
|
|
/// </summary>
|
|||
|
|
/// <param name="period">Rolling window size for the test (must be ≥ 20).</param>
|
|||
|
|
/// <param name="maxLag">Maximum number of augmented lag terms. 0 = auto-select via AIC using floor(12·(n/100)^0.25).</param>
|
|||
|
|
/// <param name="regression">Regression model type: NoConstant, Constant (default), or ConstantAndTrend.</param>
|
|||
|
|
public Adf(int period = 50, int maxLag = 0, AdfRegression regression = AdfRegression.Constant)
|
|||
|
|
{
|
|||
|
|
if (period < MinPeriod)
|
|||
|
|
{
|
|||
|
|
throw new ArgumentOutOfRangeException(nameof(period),
|
|||
|
|
$"Period must be greater than or equal to {MinPeriod} for ADF test.");
|
|||
|
|
}
|
|||
|
|
if (maxLag < 0)
|
|||
|
|
{
|
|||
|
|
throw new ArgumentOutOfRangeException(nameof(maxLag),
|
|||
|
|
"maxLag must be non-negative (0 = auto-select).");
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
_period = period;
|
|||
|
|
_maxLag = maxLag;
|
|||
|
|
_regression = regression;
|
|||
|
|
_buffer = new RingBuffer(period);
|
|||
|
|
|
|||
|
|
string regStr = regression switch
|
|||
|
|
{
|
|||
|
|
AdfRegression.NoConstant => "nc",
|
|||
|
|
AdfRegression.Constant => "c",
|
|||
|
|
AdfRegression.ConstantAndTrend => "ct",
|
|||
|
|
_ => "c"
|
|||
|
|
};
|
|||
|
|
Name = $"ADF({period},{maxLag},{regStr})";
|
|||
|
|
WarmupPeriod = period;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|||
|
|
public override TValue Update(TValue input, bool isNew = true)
|
|||
|
|
{
|
|||
|
|
double value = input.Value;
|
|||
|
|
|
|||
|
|
if (!double.IsFinite(value))
|
|||
|
|
{
|
|||
|
|
value = _lastValidValue;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
_lastValidValue = value;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if (isNew)
|
|||
|
|
{
|
|||
|
|
_inputCountSaved = _inputCount;
|
|||
|
|
_buffer.Add(value);
|
|||
|
|
_inputCount++;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
_inputCount = _inputCountSaved;
|
|||
|
|
_buffer.UpdateNewest(value);
|
|||
|
|
_inputCount++;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double pValue;
|
|||
|
|
if (_inputCount < MinPeriod)
|
|||
|
|
{
|
|||
|
|
pValue = 1.0; // Not enough data — assume unit root
|
|||
|
|
Statistic = 0.0;
|
|||
|
|
PValue = 1.0;
|
|||
|
|
LagsUsed = 0;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
var span = _buffer.GetSpan();
|
|||
|
|
var result = ComputeAdf(span, _maxLag, _regression);
|
|||
|
|
Statistic = result.Statistic;
|
|||
|
|
PValue = result.PValue;
|
|||
|
|
LagsUsed = result.LagsUsed;
|
|||
|
|
pValue = result.PValue;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
Last = new TValue(input.Time, pValue);
|
|||
|
|
PubEvent(Last, isNew);
|
|||
|
|
return Last;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
public override TSeries Update(TSeries source)
|
|||
|
|
{
|
|||
|
|
if (source.Count == 0)
|
|||
|
|
{
|
|||
|
|
return [];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
int len = source.Count;
|
|||
|
|
var t = new List<long>(len);
|
|||
|
|
var v = new List<double>(len);
|
|||
|
|
CollectionsMarshal.SetCount(t, len);
|
|||
|
|
CollectionsMarshal.SetCount(v, len);
|
|||
|
|
|
|||
|
|
var tSpan = CollectionsMarshal.AsSpan(t);
|
|||
|
|
var vSpan = CollectionsMarshal.AsSpan(v);
|
|||
|
|
|
|||
|
|
Batch(source.Values, vSpan, _period, _maxLag, _regression);
|
|||
|
|
source.Times.CopyTo(tSpan);
|
|||
|
|
|
|||
|
|
// Reset running state before priming
|
|||
|
|
_buffer.Clear();
|
|||
|
|
_lastValidValue = 0;
|
|||
|
|
_inputCount = 0;
|
|||
|
|
|
|||
|
|
// Prime: replay enough bars to reconstruct internal state
|
|||
|
|
int primeStart = Math.Max(0, len - _period);
|
|||
|
|
for (int i = primeStart; i < len; i++)
|
|||
|
|
{
|
|||
|
|
Update(source[i]);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
return new TSeries(t, v);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
public override void Reset()
|
|||
|
|
{
|
|||
|
|
_buffer.Clear();
|
|||
|
|
_lastValidValue = 0;
|
|||
|
|
_inputCount = 0;
|
|||
|
|
_inputCountSaved = 0;
|
|||
|
|
Statistic = 0;
|
|||
|
|
PValue = 1.0;
|
|||
|
|
LagsUsed = 0;
|
|||
|
|
Last = default;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
public override void Prime(ReadOnlySpan<double> source, TimeSpan? step = null)
|
|||
|
|
{
|
|||
|
|
DateTime ts = DateTime.MinValue;
|
|||
|
|
foreach (double value in source)
|
|||
|
|
{
|
|||
|
|
Update(new TValue(ts, value));
|
|||
|
|
if (step.HasValue)
|
|||
|
|
{
|
|||
|
|
ts = ts.Add(step.Value);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
public static TSeries Batch(TSeries source, int period = 50, int maxLag = 0,
|
|||
|
|
AdfRegression regression = AdfRegression.Constant)
|
|||
|
|
{
|
|||
|
|
var adf = new Adf(period, maxLag, regression);
|
|||
|
|
return adf.Update(source);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|||
|
|
public static void Batch(ReadOnlySpan<double> source, Span<double> output,
|
|||
|
|
int period = 50, int maxLag = 0, AdfRegression regression = AdfRegression.Constant)
|
|||
|
|
{
|
|||
|
|
if (source.Length != output.Length)
|
|||
|
|
{
|
|||
|
|
throw new ArgumentException("Source and output must have the same length", nameof(output));
|
|||
|
|
}
|
|||
|
|
if (period < MinPeriod)
|
|||
|
|
{
|
|||
|
|
throw new ArgumentException($"Period must be >= {MinPeriod}", nameof(period));
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
int len = source.Length;
|
|||
|
|
if (len == 0)
|
|||
|
|
{
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
CalculateScalarCore(source, output, period, maxLag, regression);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
public static (TSeries Results, Adf Indicator) Calculate(TSeries source,
|
|||
|
|
int period = 50, int maxLag = 0, AdfRegression regression = AdfRegression.Constant)
|
|||
|
|
{
|
|||
|
|
var indicator = new Adf(period, maxLag, regression);
|
|||
|
|
TSeries results = indicator.Update(source);
|
|||
|
|
return (results, indicator);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// ═══════════════════════════════════════════════════════════════
|
|||
|
|
// Private implementation
|
|||
|
|
// ═══════════════════════════════════════════════════════════════
|
|||
|
|
|
|||
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|||
|
|
private static void CalculateScalarCore(ReadOnlySpan<double> source, Span<double> output,
|
|||
|
|
int period, int maxLag, AdfRegression regression)
|
|||
|
|
{
|
|||
|
|
int len = source.Length;
|
|||
|
|
|
|||
|
|
// We need a sliding window of 'period' values
|
|||
|
|
const int StackallocThreshold = 256;
|
|||
|
|
double[]? rentedBuf = null;
|
|||
|
|
scoped Span<double> windowBuf;
|
|||
|
|
if (period <= StackallocThreshold)
|
|||
|
|
{
|
|||
|
|
windowBuf = stackalloc double[period];
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
rentedBuf = ArrayPool<double>.Shared.Rent(period);
|
|||
|
|
windowBuf = rentedBuf.AsSpan(0, period);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
try
|
|||
|
|
{
|
|||
|
|
for (int i = 0; i < len; i++)
|
|||
|
|
{
|
|||
|
|
if (i < period - 1)
|
|||
|
|
{
|
|||
|
|
output[i] = 1.0; // Not enough data
|
|||
|
|
continue;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Fill window
|
|||
|
|
int windowStart = i - period + 1;
|
|||
|
|
double pv = 0;
|
|||
|
|
for (int j = 0; j < period; j++)
|
|||
|
|
{
|
|||
|
|
double v = source[windowStart + j];
|
|||
|
|
if (!double.IsFinite(v))
|
|||
|
|
{
|
|||
|
|
v = pv;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
pv = v;
|
|||
|
|
}
|
|||
|
|
windowBuf[j] = v;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
var result = ComputeAdf(windowBuf[..period], maxLag, regression);
|
|||
|
|
output[i] = result.PValue;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
finally
|
|||
|
|
{
|
|||
|
|
if (rentedBuf is not null)
|
|||
|
|
{
|
|||
|
|
ArrayPool<double>.Shared.Return(rentedBuf);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/// <summary>
|
|||
|
|
/// Core ADF computation on a window of prices. Returns (Statistic, PValue, LagsUsed).
|
|||
|
|
/// </summary>
|
|||
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|||
|
|
private static (double Statistic, double PValue, int LagsUsed) ComputeAdf(
|
|||
|
|
ReadOnlySpan<double> y, int maxLag, AdfRegression regression)
|
|||
|
|
{
|
|||
|
|
int n = y.Length;
|
|||
|
|
if (n < 3)
|
|||
|
|
{
|
|||
|
|
return (0, 1.0, 0);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Compute first differences: Δy_t = y_t - y_{t-1}
|
|||
|
|
int nDiff = n - 1;
|
|||
|
|
const int StackallocThreshold = 256;
|
|||
|
|
|
|||
|
|
double[]? rentedDy = null;
|
|||
|
|
scoped Span<double> dy;
|
|||
|
|
if (nDiff <= StackallocThreshold)
|
|||
|
|
{
|
|||
|
|
dy = stackalloc double[nDiff];
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
rentedDy = ArrayPool<double>.Shared.Rent(nDiff);
|
|||
|
|
dy = rentedDy.AsSpan(0, nDiff);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
try
|
|||
|
|
{
|
|||
|
|
for (int i = 0; i < nDiff; i++)
|
|||
|
|
{
|
|||
|
|
dy[i] = y[i + 1] - y[i];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Auto-lag selection: Schwert (1989) rule
|
|||
|
|
int autoMaxLag = maxLag > 0
|
|||
|
|
? maxLag
|
|||
|
|
: Math.Max(1, (int)Math.Floor(12.0 * Math.Pow(n / 100.0, 0.25)));
|
|||
|
|
|
|||
|
|
// Cap lag to avoid underdetermined system
|
|||
|
|
int extraRegressors = regression switch
|
|||
|
|
{
|
|||
|
|
AdfRegression.NoConstant => 1, // γ only
|
|||
|
|
AdfRegression.Constant => 2, // α, γ
|
|||
|
|
AdfRegression.ConstantAndTrend => 3, // α, β, γ
|
|||
|
|
_ => 2
|
|||
|
|
};
|
|||
|
|
|
|||
|
|
// We need at least k+1 observations for k regressors
|
|||
|
|
int maxFeasibleLag = Math.Max(0, nDiff - extraRegressors - 2);
|
|||
|
|
autoMaxLag = Math.Min(autoMaxLag, maxFeasibleLag);
|
|||
|
|
autoMaxLag = Math.Min(autoMaxLag, MaxRegressors - extraRegressors);
|
|||
|
|
|
|||
|
|
// If auto-selecting, try each lag and pick best AIC
|
|||
|
|
int bestLag;
|
|||
|
|
if (maxLag == 0 && autoMaxLag > 0)
|
|||
|
|
{
|
|||
|
|
bestLag = SelectLagByAic(y, dy, autoMaxLag, regression);
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
bestLag = autoMaxLag;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Run ADF regression with selected lag
|
|||
|
|
var (tStat, pVal) = RunAdfRegression(y, dy, bestLag, regression, n);
|
|||
|
|
return (tStat, pVal, bestLag);
|
|||
|
|
}
|
|||
|
|
finally
|
|||
|
|
{
|
|||
|
|
if (rentedDy is not null)
|
|||
|
|
{
|
|||
|
|
ArrayPool<double>.Shared.Return(rentedDy);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/// <summary>
|
|||
|
|
/// Select optimal lag via AIC: AIC = n·ln(RSS/n) + 2·k
|
|||
|
|
/// </summary>
|
|||
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|||
|
|
private static int SelectLagByAic(ReadOnlySpan<double> y, ReadOnlySpan<double> dy,
|
|||
|
|
int maxLag, AdfRegression regression)
|
|||
|
|
{
|
|||
|
|
double bestAic = double.MaxValue;
|
|||
|
|
int bestLag = 0;
|
|||
|
|
int nOrig = y.Length;
|
|||
|
|
|
|||
|
|
for (int lag = 0; lag <= maxLag; lag++)
|
|||
|
|
{
|
|||
|
|
int nObs = dy.Length - lag;
|
|||
|
|
if (nObs < 3)
|
|||
|
|
{
|
|||
|
|
break;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
int k = (regression switch
|
|||
|
|
{
|
|||
|
|
AdfRegression.NoConstant => 1,
|
|||
|
|
AdfRegression.Constant => 2,
|
|||
|
|
AdfRegression.ConstantAndTrend => 3,
|
|||
|
|
_ => 2
|
|||
|
|
}) + lag;
|
|||
|
|
|
|||
|
|
if (nObs <= k + 1)
|
|||
|
|
{
|
|||
|
|
break;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double rss = ComputeRss(y, dy, lag, regression, nOrig);
|
|||
|
|
|
|||
|
|
if (rss <= 0 || !double.IsFinite(rss))
|
|||
|
|
{
|
|||
|
|
continue;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double aic = nObs * Math.Log(rss / nObs) + 2.0 * k;
|
|||
|
|
if (aic < bestAic)
|
|||
|
|
{
|
|||
|
|
bestAic = aic;
|
|||
|
|
bestLag = lag;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
return bestLag;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/// <summary>
|
|||
|
|
/// Compute residual sum of squares for a given lag configuration.
|
|||
|
|
/// </summary>
|
|||
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|||
|
|
private static double ComputeRss(ReadOnlySpan<double> y, ReadOnlySpan<double> dy,
|
|||
|
|
int lag, AdfRegression regression, int nOrig)
|
|||
|
|
{
|
|||
|
|
int nDiff = dy.Length;
|
|||
|
|
int startIdx = lag;
|
|||
|
|
int nObs = nDiff - startIdx;
|
|||
|
|
if (nObs < 3)
|
|||
|
|
{
|
|||
|
|
return double.MaxValue;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
int k = (regression switch
|
|||
|
|
{
|
|||
|
|
AdfRegression.NoConstant => 1,
|
|||
|
|
AdfRegression.Constant => 2,
|
|||
|
|
AdfRegression.ConstantAndTrend => 3,
|
|||
|
|
_ => 2
|
|||
|
|
}) + lag;
|
|||
|
|
|
|||
|
|
if (k > MaxRegressors || nObs <= k)
|
|||
|
|
{
|
|||
|
|
return double.MaxValue;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Build X'X and X'y via accumulation (no explicit matrix allocation)
|
|||
|
|
Span<double> xtx = stackalloc double[k * k];
|
|||
|
|
Span<double> xty = stackalloc double[k];
|
|||
|
|
xtx.Clear();
|
|||
|
|
xty.Clear();
|
|||
|
|
|
|||
|
|
Span<double> xRow = stackalloc double[k];
|
|||
|
|
|
|||
|
|
for (int t = startIdx; t < nDiff; t++)
|
|||
|
|
{
|
|||
|
|
// Build regressor row
|
|||
|
|
int col = 0;
|
|||
|
|
|
|||
|
|
// Deterministic terms
|
|||
|
|
if (regression is AdfRegression.Constant or AdfRegression.ConstantAndTrend)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = 1.0; // intercept
|
|||
|
|
}
|
|||
|
|
if (regression == AdfRegression.ConstantAndTrend)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = t + 1; // trend
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// y_{t-1} (the lagged level — this is the key regressor)
|
|||
|
|
xRow[col++] = y[t]; // y[t] in 0-based corresponds to y_{t} which is lagged level
|
|||
|
|
|
|||
|
|
// Lagged differences
|
|||
|
|
for (int j = 1; j <= lag; j++)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = dy[t - j];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double yVal = dy[t];
|
|||
|
|
|
|||
|
|
// Accumulate X'X and X'y
|
|||
|
|
for (int r = 0; r < k; r++)
|
|||
|
|
{
|
|||
|
|
for (int c = r; c < k; c++)
|
|||
|
|
{
|
|||
|
|
xtx[r * k + c] += xRow[r] * xRow[c];
|
|||
|
|
}
|
|||
|
|
xty[r] += xRow[r] * yVal;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Fill lower triangle of symmetric X'X
|
|||
|
|
for (int r = 1; r < k; r++)
|
|||
|
|
{
|
|||
|
|
for (int c = 0; c < r; c++)
|
|||
|
|
{
|
|||
|
|
xtx[r * k + c] = xtx[c * k + r];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Solve via Cholesky decomposition
|
|||
|
|
Span<double> beta = stackalloc double[k];
|
|||
|
|
if (!SolveCholesky(xtx, xty, beta, k))
|
|||
|
|
{
|
|||
|
|
return double.MaxValue;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Compute RSS
|
|||
|
|
double rss = 0;
|
|||
|
|
for (int t = startIdx; t < nDiff; t++)
|
|||
|
|
{
|
|||
|
|
int col = 0;
|
|||
|
|
if (regression is AdfRegression.Constant or AdfRegression.ConstantAndTrend)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = 1.0;
|
|||
|
|
}
|
|||
|
|
if (regression == AdfRegression.ConstantAndTrend)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = t + 1;
|
|||
|
|
}
|
|||
|
|
xRow[col++] = y[t];
|
|||
|
|
for (int j = 1; j <= lag; j++)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = dy[t - j];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double predicted = 0;
|
|||
|
|
for (int c = 0; c < k; c++)
|
|||
|
|
{
|
|||
|
|
predicted += beta[c] * xRow[c];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double residual = dy[t] - predicted;
|
|||
|
|
rss += residual * residual;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
return rss;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/// <summary>
|
|||
|
|
/// Run ADF regression and return (t-statistic, p-value).
|
|||
|
|
/// </summary>
|
|||
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|||
|
|
private static (double tStat, double pValue) RunAdfRegression(
|
|||
|
|
ReadOnlySpan<double> y, ReadOnlySpan<double> dy,
|
|||
|
|
int lag, AdfRegression regression, int nOrig)
|
|||
|
|
{
|
|||
|
|
int nDiff = dy.Length;
|
|||
|
|
int startIdx = lag;
|
|||
|
|
int nObs = nDiff - startIdx;
|
|||
|
|
if (nObs < 3)
|
|||
|
|
{
|
|||
|
|
return (0, 1.0);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
int k = (regression switch
|
|||
|
|
{
|
|||
|
|
AdfRegression.NoConstant => 1,
|
|||
|
|
AdfRegression.Constant => 2,
|
|||
|
|
AdfRegression.ConstantAndTrend => 3,
|
|||
|
|
_ => 2
|
|||
|
|
}) + lag;
|
|||
|
|
|
|||
|
|
if (k > MaxRegressors || nObs <= k)
|
|||
|
|
{
|
|||
|
|
return (0, 1.0);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Index of γ (coefficient on y_{t-1})
|
|||
|
|
int gammaIdx = regression switch
|
|||
|
|
{
|
|||
|
|
AdfRegression.NoConstant => 0,
|
|||
|
|
AdfRegression.Constant => 1,
|
|||
|
|
AdfRegression.ConstantAndTrend => 2,
|
|||
|
|
_ => 1
|
|||
|
|
};
|
|||
|
|
|
|||
|
|
// Build X'X and X'y
|
|||
|
|
Span<double> xtx = stackalloc double[k * k];
|
|||
|
|
Span<double> xty = stackalloc double[k];
|
|||
|
|
xtx.Clear();
|
|||
|
|
xty.Clear();
|
|||
|
|
|
|||
|
|
Span<double> xRow = stackalloc double[k];
|
|||
|
|
|
|||
|
|
for (int t = startIdx; t < nDiff; t++)
|
|||
|
|
{
|
|||
|
|
int col = 0;
|
|||
|
|
if (regression is AdfRegression.Constant or AdfRegression.ConstantAndTrend)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = 1.0;
|
|||
|
|
}
|
|||
|
|
if (regression == AdfRegression.ConstantAndTrend)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = t + 1;
|
|||
|
|
}
|
|||
|
|
xRow[col++] = y[t];
|
|||
|
|
for (int j = 1; j <= lag; j++)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = dy[t - j];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double yVal = dy[t];
|
|||
|
|
|
|||
|
|
for (int r = 0; r < k; r++)
|
|||
|
|
{
|
|||
|
|
for (int c = r; c < k; c++)
|
|||
|
|
{
|
|||
|
|
xtx[r * k + c] += xRow[r] * xRow[c];
|
|||
|
|
}
|
|||
|
|
xty[r] += xRow[r] * yVal;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Fill lower triangle
|
|||
|
|
for (int r = 1; r < k; r++)
|
|||
|
|
{
|
|||
|
|
for (int c = 0; c < r; c++)
|
|||
|
|
{
|
|||
|
|
xtx[r * k + c] = xtx[c * k + r];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Solve for beta via Cholesky
|
|||
|
|
Span<double> beta = stackalloc double[k];
|
|||
|
|
if (!SolveCholesky(xtx, xty, beta, k))
|
|||
|
|
{
|
|||
|
|
return (0, 1.0);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Compute (X'X)^-1 for SE calculation — rebuild and invert
|
|||
|
|
xtx.Clear();
|
|||
|
|
xty.Clear();
|
|||
|
|
|
|||
|
|
// Rebuild X'X
|
|||
|
|
for (int t = startIdx; t < nDiff; t++)
|
|||
|
|
{
|
|||
|
|
int col = 0;
|
|||
|
|
if (regression is AdfRegression.Constant or AdfRegression.ConstantAndTrend)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = 1.0;
|
|||
|
|
}
|
|||
|
|
if (regression == AdfRegression.ConstantAndTrend)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = t + 1;
|
|||
|
|
}
|
|||
|
|
xRow[col++] = y[t];
|
|||
|
|
for (int j = 1; j <= lag; j++)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = dy[t - j];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
for (int r = 0; r < k; r++)
|
|||
|
|
{
|
|||
|
|
for (int c = r; c < k; c++)
|
|||
|
|
{
|
|||
|
|
xtx[r * k + c] += xRow[r] * xRow[c];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
for (int r = 1; r < k; r++)
|
|||
|
|
{
|
|||
|
|
for (int c = 0; c < r; c++)
|
|||
|
|
{
|
|||
|
|
xtx[r * k + c] = xtx[c * k + r];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Compute RSS
|
|||
|
|
double rss = 0;
|
|||
|
|
for (int t = startIdx; t < nDiff; t++)
|
|||
|
|
{
|
|||
|
|
int col = 0;
|
|||
|
|
if (regression is AdfRegression.Constant or AdfRegression.ConstantAndTrend)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = 1.0;
|
|||
|
|
}
|
|||
|
|
if (regression == AdfRegression.ConstantAndTrend)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = t + 1;
|
|||
|
|
}
|
|||
|
|
xRow[col++] = y[t];
|
|||
|
|
for (int j = 1; j <= lag; j++)
|
|||
|
|
{
|
|||
|
|
xRow[col++] = dy[t - j];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double predicted = 0;
|
|||
|
|
for (int c = 0; c < k; c++)
|
|||
|
|
{
|
|||
|
|
predicted += beta[c] * xRow[c];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double residual = dy[t] - predicted;
|
|||
|
|
rss += residual * residual;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double s2 = rss / (nObs - k); // residual variance estimate
|
|||
|
|
if (s2 <= 0 || !double.IsFinite(s2))
|
|||
|
|
{
|
|||
|
|
return (0, 1.0);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Invert X'X via Cholesky to get (X'X)^-1
|
|||
|
|
Span<double> xtxInv = stackalloc double[k * k];
|
|||
|
|
if (!InvertViaCholesky(xtx, xtxInv, k))
|
|||
|
|
{
|
|||
|
|
return (0, 1.0);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double seGamma = Math.Sqrt(s2 * xtxInv[gammaIdx * k + gammaIdx]);
|
|||
|
|
if (seGamma <= 1e-15 || !double.IsFinite(seGamma))
|
|||
|
|
{
|
|||
|
|
return (0, 1.0);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double tStat = beta[gammaIdx] / seGamma;
|
|||
|
|
double pVal = MacKinnonPValue(tStat, regression, nOrig);
|
|||
|
|
|
|||
|
|
return (tStat, pVal);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// ═══════════════════════════════════════════════════════════════
|
|||
|
|
// Cholesky solver for k×k system (k ≤ 8, stackalloc safe)
|
|||
|
|
// ═══════════════════════════════════════════════════════════════
|
|||
|
|
|
|||
|
|
/// <summary>
|
|||
|
|
/// Solves A·x = b via Cholesky decomposition where A is symmetric positive definite.
|
|||
|
|
/// Returns false if decomposition fails (matrix not positive definite).
|
|||
|
|
/// </summary>
|
|||
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|||
|
|
private static bool SolveCholesky(ReadOnlySpan<double> A, ReadOnlySpan<double> b,
|
|||
|
|
Span<double> x, int k)
|
|||
|
|
{
|
|||
|
|
// Cholesky decomposition: A = L·L'
|
|||
|
|
Span<double> L = stackalloc double[k * k];
|
|||
|
|
L.Clear();
|
|||
|
|
|
|||
|
|
for (int i = 0; i < k; i++)
|
|||
|
|
{
|
|||
|
|
for (int j = 0; j <= i; j++)
|
|||
|
|
{
|
|||
|
|
double sum = A[i * k + j];
|
|||
|
|
for (int p = 0; p < j; p++)
|
|||
|
|
{
|
|||
|
|
sum -= L[i * k + p] * L[j * k + p];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if (i == j)
|
|||
|
|
{
|
|||
|
|
if (sum <= 1e-15)
|
|||
|
|
{
|
|||
|
|
return false; // Not positive definite
|
|||
|
|
}
|
|||
|
|
L[i * k + j] = Math.Sqrt(sum);
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
L[i * k + j] = sum / L[j * k + j];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Forward substitution: L·z = b
|
|||
|
|
Span<double> z = stackalloc double[k];
|
|||
|
|
for (int i = 0; i < k; i++)
|
|||
|
|
{
|
|||
|
|
double sum = b[i];
|
|||
|
|
for (int j = 0; j < i; j++)
|
|||
|
|
{
|
|||
|
|
sum -= L[i * k + j] * z[j];
|
|||
|
|
}
|
|||
|
|
z[i] = sum / L[i * k + i];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Back substitution: L'·x = z
|
|||
|
|
for (int i = k - 1; i >= 0; i--)
|
|||
|
|
{
|
|||
|
|
double sum = z[i];
|
|||
|
|
for (int j = i + 1; j < k; j++)
|
|||
|
|
{
|
|||
|
|
sum -= L[j * k + i] * x[j];
|
|||
|
|
}
|
|||
|
|
x[i] = sum / L[i * k + i];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
return true;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/// <summary>
|
|||
|
|
/// Inverts a symmetric positive definite matrix via Cholesky decomposition.
|
|||
|
|
/// </summary>
|
|||
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|||
|
|
private static bool InvertViaCholesky(ReadOnlySpan<double> A, Span<double> Ainv, int k)
|
|||
|
|
{
|
|||
|
|
// Cholesky: A = L·L'
|
|||
|
|
Span<double> L = stackalloc double[k * k];
|
|||
|
|
L.Clear();
|
|||
|
|
|
|||
|
|
for (int i = 0; i < k; i++)
|
|||
|
|
{
|
|||
|
|
for (int j = 0; j <= i; j++)
|
|||
|
|
{
|
|||
|
|
double sum = A[i * k + j];
|
|||
|
|
for (int p = 0; p < j; p++)
|
|||
|
|
{
|
|||
|
|
sum -= L[i * k + p] * L[j * k + p];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if (i == j)
|
|||
|
|
{
|
|||
|
|
if (sum <= 1e-15)
|
|||
|
|
{
|
|||
|
|
return false;
|
|||
|
|
}
|
|||
|
|
L[i * k + j] = Math.Sqrt(sum);
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
L[i * k + j] = sum / L[j * k + j];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Invert L (lower triangular)
|
|||
|
|
Span<double> Linv = stackalloc double[k * k];
|
|||
|
|
Linv.Clear();
|
|||
|
|
for (int i = 0; i < k; i++)
|
|||
|
|
{
|
|||
|
|
Linv[i * k + i] = 1.0 / L[i * k + i];
|
|||
|
|
for (int j = i + 1; j < k; j++)
|
|||
|
|
{
|
|||
|
|
double sum = 0;
|
|||
|
|
for (int p = i; p < j; p++)
|
|||
|
|
{
|
|||
|
|
sum -= L[j * k + p] * Linv[p * k + i];
|
|||
|
|
}
|
|||
|
|
Linv[j * k + i] = sum / L[j * k + j];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// A^-1 = (L')^-1 · L^-1 = Linv' · Linv
|
|||
|
|
Ainv.Clear();
|
|||
|
|
for (int i = 0; i < k; i++)
|
|||
|
|
{
|
|||
|
|
for (int j = i; j < k; j++)
|
|||
|
|
{
|
|||
|
|
double sum = 0;
|
|||
|
|
for (int p = j; p < k; p++)
|
|||
|
|
{
|
|||
|
|
sum += Linv[p * k + i] * Linv[p * k + j];
|
|||
|
|
}
|
|||
|
|
Ainv[i * k + j] = sum;
|
|||
|
|
Ainv[j * k + i] = sum;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
return true;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// ═══════════════════════════════════════════════════════════════
|
|||
|
|
// MacKinnon (1994) p-value approximation
|
|||
|
|
// ═══════════════════════════════════════════════════════════════
|
|||
|
|
|
|||
|
|
// MacKinnon (2010) coefficients for N=1 (univariate)
|
|||
|
|
// [tau_star, tau_min, tau_max]
|
|||
|
|
// SmallP: polynomial coefficients for small-p region (degree 2)
|
|||
|
|
// LargeP: polynomial coefficients for large-p region (degree 3)
|
|||
|
|
|
|||
|
|
// NoConstant
|
|||
|
|
private static readonly double[] SmallP_Nc = [0.6344, 1.2378, 0.032496];
|
|||
|
|
private static readonly double[] LargeP_Nc = [0.4797, 0.93557, -0.06999, 0.033066];
|
|||
|
|
|
|||
|
|
// Constant
|
|||
|
|
private static readonly double[] SmallP_C = [2.1659, 1.4412, 0.038269];
|
|||
|
|
private static readonly double[] LargeP_C = [1.7339, 0.93202, -0.12745, -0.010368];
|
|||
|
|
|
|||
|
|
// ConstantAndTrend
|
|||
|
|
private static readonly double[] SmallP_Ct = [3.2512, 1.6047, 0.049588];
|
|||
|
|
private static readonly double[] LargeP_Ct = [2.5261, 0.61654, -0.37956, -0.060285];
|
|||
|
|
|
|||
|
|
/// <summary>
|
|||
|
|
/// Computes the MacKinnon approximate p-value for the ADF test statistic.
|
|||
|
|
/// Uses interpolation coefficients from MacKinnon (1994, 2010) for N=1.
|
|||
|
|
/// </summary>
|
|||
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|||
|
|
private static double MacKinnonPValue(double tStat, AdfRegression regression, int n)
|
|||
|
|
{
|
|||
|
|
double tauStar, tauMin, tauMax;
|
|||
|
|
double[] smallP, largeP;
|
|||
|
|
|
|||
|
|
switch (regression)
|
|||
|
|
{
|
|||
|
|
case AdfRegression.NoConstant:
|
|||
|
|
tauStar = -1.04;
|
|||
|
|
tauMin = -19.04;
|
|||
|
|
tauMax = double.PositiveInfinity;
|
|||
|
|
smallP = SmallP_Nc;
|
|||
|
|
largeP = LargeP_Nc;
|
|||
|
|
break;
|
|||
|
|
|
|||
|
|
case AdfRegression.ConstantAndTrend:
|
|||
|
|
tauStar = -2.89;
|
|||
|
|
tauMin = -16.18;
|
|||
|
|
tauMax = 0.70;
|
|||
|
|
smallP = SmallP_Ct;
|
|||
|
|
largeP = LargeP_Ct;
|
|||
|
|
break;
|
|||
|
|
|
|||
|
|
default: // Constant
|
|||
|
|
tauStar = -1.61;
|
|||
|
|
tauMin = -18.83;
|
|||
|
|
tauMax = 2.74;
|
|||
|
|
smallP = SmallP_C;
|
|||
|
|
largeP = LargeP_C;
|
|||
|
|
break;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Clamp to valid range
|
|||
|
|
if (tStat < tauMin)
|
|||
|
|
{
|
|||
|
|
return 0.0; // Extremely stationary
|
|||
|
|
}
|
|||
|
|
if (double.IsFinite(tauMax) && tStat > tauMax)
|
|||
|
|
{
|
|||
|
|
return 1.0;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if (tStat <= tauStar)
|
|||
|
|
{
|
|||
|
|
// Small p-value region: p = NormCdf(poly3(tStat))
|
|||
|
|
double poly = smallP[0] + smallP[1] * tStat + smallP[2] * tStat * tStat;
|
|||
|
|
return NormCdf(poly);
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
// Large p-value region: p = NormCdf(poly4(tStat))
|
|||
|
|
double t2 = tStat * tStat;
|
|||
|
|
double t3 = t2 * tStat;
|
|||
|
|
double poly = largeP[0] + largeP[1] * tStat + largeP[2] * t2 + largeP[3] * t3;
|
|||
|
|
return NormCdf(poly);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/// <summary>
|
|||
|
|
/// Standard normal CDF: Φ(x) = 0.5 * erfc(-x / √2)
|
|||
|
|
/// Machine precision via .NET BCL Math.Erfc.
|
|||
|
|
/// </summary>
|
|||
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|||
|
|
private static double NormCdf(double x)
|
|||
|
|
{
|
|||
|
|
// erfc is not in Math directly; use the identity:
|
|||
|
|
// Φ(x) = 0.5 * (1 + erf(x / √2))
|
|||
|
|
// But .NET doesn't have erf either. Use the rational approximation instead.
|
|||
|
|
return NormCdfApprox(x);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/// <summary>
|
|||
|
|
/// Abramowitz and Stegun 7.1.26 erf approximation → standard normal CDF.
|
|||
|
|
/// Φ(x) = 0.5·(1 + erf(x/√2)), where erf is computed via the A&S rational
|
|||
|
|
/// approximation with coefficients a1…a5, p = 0.3275911.
|
|||
|
|
/// Maximum error: 7.5e-8.
|
|||
|
|
/// </summary>
|
|||
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|||
|
|
private static double NormCdfApprox(double x)
|
|||
|
|
{
|
|||
|
|
const double a1 = 0.254829592;
|
|||
|
|
const double a2 = -0.284496736;
|
|||
|
|
const double a3 = 1.421413741;
|
|||
|
|
const double a4 = -1.453152027;
|
|||
|
|
const double a5 = 1.061405429;
|
|||
|
|
const double p = 0.3275911;
|
|||
|
|
const double inv_sqrt2 = 0.70710678118654752;
|
|||
|
|
|
|||
|
|
int sign = x < 0 ? -1 : 1;
|
|||
|
|
|
|||
|
|
// erf argument: z = |x| / √2
|
|||
|
|
double z = Math.Abs(x) * inv_sqrt2;
|
|||
|
|
|
|||
|
|
// Rational approximation for erfc(z) = poly · exp(-z²)
|
|||
|
|
double t = 1.0 / Math.FusedMultiplyAdd(p, z, 1.0);
|
|||
|
|
double t2 = t * t;
|
|||
|
|
double t3 = t2 * t;
|
|||
|
|
double t4 = t3 * t;
|
|||
|
|
double t5 = t4 * t;
|
|||
|
|
|
|||
|
|
double poly = a1 * t + a2 * t2 + a3 * t3 + a4 * t4 + a5 * t5;
|
|||
|
|
double erfcZ = poly * Math.Exp(-z * z);
|
|||
|
|
double erfZ = 1.0 - erfcZ;
|
|||
|
|
|
|||
|
|
// Φ(x) = 0.5 · (1 + sign · erf(|x|/√2))
|
|||
|
|
return 0.5 * (1.0 + sign * erfZ);
|
|||
|
|
}
|
|||
|
|
}
|