WIPFunctionLyaponov

// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © RicardoSantos

//@version=5
// checkpoint
// 100 matrix version
// 126 map version

//@description Lyapunov exponents are mathematical measures used to describe the behavior of a system over
// time. They are named after Russian mathematician Alexei Lyapunov, who first introduced the concept in the
// late 19th century. The exponent is defined as the rate at which a particular function or variable changes
// over time, and can be positive, negative, or zero.
//  Positive exponents indicate that a system tends to grow or expand over time, while negative exponents
// indicate that a system tends to shrink or decay. Zero exponents indicate that the system does not change
// significantly over time. Lyapunov exponents are used in various fields of science and engineering, including
// physics, economics, and biology, to study the long-term behavior of complex systems.
// ~ generated description from vicuna13b
//
// ---
// To calculate the Lyapunov Exponent (LE) of a given Time Series, we need to follow these steps:
// 1. Firstly, you should have access to your data in some format like CSV or Excel file. If not, then you can collect it manually using tools such as stopwatches and measuring tapes.
// 2. Once the data is collected, clean it up by removing any outliers that may skew results. This step involves checking for inconsistencies within your dataset (e.g., extremely large or small values) and either discarding them entirely or replacing with more reasonable estimates based on surrounding values.
// 3. Next, you need to determine the dimension of your time series data. In most cases, this will be equal to the number of variables being measured in each observation period (e.g., temperature, humidity, wind speed).
// 4. Now that we have a clean dataset with known dimensions, we can calculate the LE for our Time Series using the following formula:
//    λ = log(||M^T * M - I||)/log(||v||)
//    where:
//      λ (Lyapunov Exponent) is the quantity that will be calculated.
//      ||...|| denotes an Euclidean norm of a vector or matrix, which essentially means taking the square root of the sum of squares for each element in the vector/matrix.
//     M represents our Jacobian Matrix whose elements are given by:
//        J_ij = (∂fj / ∂xj) where fj is the jth variable and xj is the ith component of the initial condition vector x(t). In other words, each element in this matrix represents how much a small change in one variable affects another.
//     I denotes an identity matrix whose elements are all equal to 1 (or any constant value if you prefer). This term essentially acts as a baseline for comparison purposes since we want our Jacobian Matrix M^T * M to be close to it when the system is stable and far away from it when the system is unstable.
//     v represents an arbitrary vector whose Euclidean norm ||v|| will serve as a scaling factor in our calculation. The choice of this particular vector does not matter since we are only interested in its magnitude (i.e., length) for purposes of normalization. However, if you want to ensure that your results are accurate and consistent across different datasets or scenarios, it is recommended to use the same initial condition vector x(t) as used earlier when calculating our Jacobian Matrix M.
// 5. Finally, once we have calculated λ using the formula above, we can interpret its value in terms of stability/instability for our Time Series data:
//     - If λ < 0, then this indicates that the system is stable (i.e., nearby trajectories will converge towards each other over time).
//     - On the other hand, if λ > 0, then this implies that the system is unstable (i.e., nearby trajectories will diverge away from one another over time).
// ~ generated description from airoboros33b
//
// ---
// Reference:
// https://en.wikipedia.org/wiki/Lyapunov_exponent
// https://www.collimator.ai/reference-guides/what-is-a-lyapunov-exponent
// https://blog.abhranil.net/2014/07/22/calculating-the-lyapunov-exponent-of-a-time-series-with-python-code/
// https://www.researchgate.net/publication/222460963_Determining_Lyapunov_Exponents_From_a_Time_Series
// https://physics.stackexchange.com/questions/102529/calculating-lyapunov-exponents-from-a-multi-dimensional-experimental-time-series
//
// ---
// This is a work in progress, it may contain errors so use with caution.
// If you find flaws or suggest something new, please leave a comment bellow.
library("WIPFunctionLyaponov")

import RicardoSantos/SimilarityMeasures/1
import RicardoSantos/CommonTypesMapUtil/1

// @function helper function to get the name of distance function by a index (0 -> 13).\
// Functions: SSD, Euclidean, Manhattan, Minkowski, Chebyshev, Correlation, Cosine, Camberra, MAE, MSE, Lorentzian, Intersection, Penrose Shape, Meehl.
export _measure_function (int i=0) =>
	switch i
		00 => 'ssd'
		01 => 'euclidean'
		02 => 'manhattan'
		03 => 'minkowski'
		04 => 'chebychev'
		05 => 'correlation'
		06 => 'cosine'
		07 => 'camberra'
		08 => 'mae'
		09 => 'mse'
		10 => 'lorentzian'
		11 => 'intersection'
		12 => 'penrose'
		13 => 'meehl'

// @function Calculate the vector distance between two vectors.
vdistance (array<float> a, array<float> b, string type = 'ssd') =>
	switch type
		'ssd'       	=> SimilarityMeasures.ssd(a, b)
		'euclidean' 	=> SimilarityMeasures.euclidean(a, b)
		'manhattan' 	=> SimilarityMeasures.manhattan(a, b)
		'minkowski' 	=> SimilarityMeasures.minkowski(a, b)
		'chebychev' 	=> SimilarityMeasures.chebyshev(a, b)
		'correlation' 	=> SimilarityMeasures.correlation(a, b)
		'cosine' 		=> SimilarityMeasures.cosine(a, b)
		'camberra' 		=> SimilarityMeasures.camberra(a, b)
		'mae' 			=> SimilarityMeasures.mae(a, b)
		'mse' 			=> SimilarityMeasures.mse(a, b)
		'lorentzian' 	=> SimilarityMeasures.lorentzian(a, b)
		'intersection' 	=> SimilarityMeasures.intersection(a, b)
		'penrose' 		=> SimilarityMeasures.penrose(a, b)
		'meehl' 		=> SimilarityMeasures.meehl(a, b)

// @function Helper function to test the output exponents state system and outputs description into a string.
export _test (array<float> L) =>
	float _sum = 0.0
	int _neg_count = 0
	bool _is_chaotic = false
	for _e in L
		if _e > 0.0
			_is_chaotic := true
			break
		if _e < 0.0
			_neg_count += 1
		_sum += _e
	_is_stable = _neg_count == L.size()
	switch
		_is_chaotic => 'System is Chaotic.\n Small differences in initial conditions will lead to vastly different outcomes over time.\n' // atleast one positive exponent.
		_is_stable  => 'System is Stable.\n Small perturbations from equilibrium will decay over time, and the system will return to its original state.\n' // all exponents are negative
		_sum == 0.0 => 'System is Neutrally Stable.\n Small perturbations from equilibrium will neither converge nor diverge.\n'//the sum of all Lyapunov exponents is zero
		=> 'Unable to determine System.'

// reference: arXiv:2308.01013
// Bayesian framework for characterizing
// cryptocurrency market dynamics, structural
// dependency, and volatility using potential field
// by: Anoop C V, Neeraj Negi, Anup Aprem


// @function Estimate the Lyaponov Exponents for multiple series in a row matrix.
// @param data				Data series in a row matrix format (a row represents a 1d series).
// @param initial_distance	Initial distance limit.
// @param distance_function	Name of the distance function to be used, default:`ssd`.
// @returns List of Lyaponov exponents.
export estimate (map<string, CommonTypesMapUtil.ArrayFloat> X, float initial_distance=1.0e6, string distance_function='ssd') =>
	array<string> _series = X.keys()
	int _N = _series.size()-1
	int _length = X.get(_series.get(0)).v.size()
	array<float> _L = array.new<float>()
	matrix<float> _lambda = matrix.new<float>(_N+1, _length, 0.0)
	for [_i, _si] in _series
		if _i < _N
			for _j = _i + 1 to _N
				array<float> _Xi = X.get(_si            ).v
				array<float> _Xj = X.get(_series.get(_j)).v
				float _dist1 = vdistance(_Xi, _Xj, distance_function)
				int _Nij = _N - _j
				// Current distance under initial distance or predetermined.
				if _dist1 <= initial_distance
					for _k = 0 to _Nij
						float _1k = 1.0 / _k
						array<float> _Sik = X.get(_series.get(_i + _k)).v
						array<float> _Sjk = X.get(_series.get(_j + _k)).v
						// Calculate divergence from initial.
						for _l = 0 to _length-1
							float _dist2 = vdistance(_Sik.slice(_l, _length), _Sjk.slice(_l, _length), distance_function)
							float _div = _1k * math.log(_dist2 / _dist1)
							_lambda.set(_k, _l, _div)
							_lambda.set(_i, _j, nz(_lambda.get(_i, _j)) + _div)
				_L.push(nz(_lambda.get(_i, _j) / _Nij, 0.0))
	_L

// @function Maximal Lyaponov Exponent.
// @param L		List of Lyapunov exponents.
// @returns Highest exponent.
export max (array<float> L) => L.max()

import RicardoSantos/DebugConsole/13 as console
logger = console.new()
logger.table.cell_set_height(0, 0, 80.0)
logger.table.cell_set_width(0, 0, 80.0)

var matrix<float> mat = matrix.new<float>(5, 10, na)
// test with random source:
// if barstate.islastconfirmedhistory
//     for _i = 0 to 4
//         for _j = 0 to 9
//             mat.set(_i, _j, 1000.0 + 100.0 * math.random(-1.0, 1.0))
//     lyap = estimate(mat, array.new<float>(10), e, k, dt)
//     log.warning('{0}: {1}', test(lyap), lyap)
//     line.new(bar_index, 0, bar_index+1, lyap.get(0))
//     for _i = 1 to 9
//         line.new(bar_index+_i, lyap.get(_i-1), bar_index+_i+1, lyap.get(_i))

// test with assets source:
string _sym0 = 'BTCUSD'
string _sym1 = 'ETHUSD'
string _sym2 = 'XRPUSD'
string _sym3 = 'SOLUSD'
string _sym4 = 'SPX'
var array<string> symbols = array.from(_sym0, _sym1, _sym2, _sym3, _sym4)
float _s0 = request.security(_sym0, timeframe.period, close, barmerge.gaps_on, barmerge.lookahead_off)
float _s1 = request.security(_sym1, timeframe.period, close, barmerge.gaps_on, barmerge.lookahead_off)
float _s2 = request.security(_sym2, timeframe.period, close, barmerge.gaps_on, barmerge.lookahead_off)
float _s3 = request.security(_sym3, timeframe.period, close, barmerge.gaps_on, barmerge.lookahead_off)
float _s4 = request.security(_sym4, timeframe.period, close, barmerge.gaps_on, barmerge.lookahead_off)

var map0 = map.new<string, CommonTypesMapUtil.ArrayFloat>()
var map1 = map.new<string, CommonTypesMapUtil.ArrayFloat>()
var map2 = map.new<string, CommonTypesMapUtil.ArrayFloat>()
// mat3 = matrix.new<float>(5, 10, 0.0)
// for _i = 0 to 4
//     for _j = 0 to 9
//         mat.set(_i, _j, nz(sel(_i, _j)))
//         mat2.set(_i, _j, bar_index[_j])
//         mat3.set(_i, _j, math.random(-1.0, 1.0))
int length = input.int(10)
if barstate.islastconfirmedhistory
	for [_si, _sym] in symbols
		_series = CommonTypesMapUtil.ArrayFloat.new(array.new<float>(length, 0.0))
		map0.put(_sym, _series)
	map1.put('seq1', CommonTypesMapUtil.ArrayFloat.new(array.new<float>(length, 0.0)))
	map1.put('seq2', CommonTypesMapUtil.ArrayFloat.new(array.new<float>(length, 0.0)))
	map1.put('seq3', CommonTypesMapUtil.ArrayFloat.new(array.new<float>(length, 0.0)))
	map2.put('rng1', CommonTypesMapUtil.ArrayFloat.new(array.new<float>(length, 0.0)))
	map2.put('rng2', CommonTypesMapUtil.ArrayFloat.new(array.new<float>(length, 0.0)))
	map2.put('rng3', CommonTypesMapUtil.ArrayFloat.new(array.new<float>(length, 0.0)))
	map2.put('rng4', CommonTypesMapUtil.ArrayFloat.new(array.new<float>(length, 0.0)))
	for _i = 0 to length - 1
		map0.get(_sym0).v.set(_i, _s0[_i])
		map0.get(_sym1).v.set(_i, _s1[_i])
		map0.get(_sym2).v.set(_i, _s2[_i])
		map0.get(_sym3).v.set(_i, _s3[_i])
		map0.get(_sym4).v.set(_i, _s4[_i])
		map1.get('seq1').v.set(_i, bar_index[_i])
		map1.get('seq2').v.set(_i, bar_index[_i])
		map1.get('seq3').v.set(_i, bar_index[_i])
		map2.get('rng1').v.set(_i, math.random(-1.0, 1.0))
		map2.get('rng2').v.set(_i, math.random(-1.5, 1.0))
		map2.get('rng3').v.set(_i, math.random(-2.0, 3.0))
		map2.get('rng4').v.set(_i, math.random(-3.0, 5.0))
	// logger.queue(str.format('{0}', map0.get(_sym0).v))
	lyap0 = estimate(map0, vdistance(map0.get(symbols.get(0)).v, map0.get(symbols.get(1)).v))
	logger.queue(str.format('{0}: {1}', _test(lyap0), str.tostring(lyap0, '#.####')))
	lyap1 = estimate(map1, vdistance(map1.get('seq1').v, map1.get('seq2').v))
	logger.queue(str.format('{0}: {1}', _test(lyap1), str.tostring(lyap1, '#.####')))
	lyap2 = estimate(map2, vdistance(map2.get('rng1').v, map2.get('rng2').v))
	logger.queue(str.format('{0}: {1}', _test(lyap2), str.tostring(lyap2, '#.####')))
//     line.new(bar_index, 0, bar_index+1, lyap.get(0))
//     for _i = 1 to 9
//         line.new(bar_index+_i, lyap.get(_i-1), bar_index+_i+1, lyap.get(_i))

logger.update()