# A Decadal Analysis of the Lead-Lag Effect in the NYSE

##  1 Introduction

As is widely known, the stock market is a complex system in which a multitude
of factors influence the performance of individual stocks and the market as a
whole. One method for comprehending—and potentially predicting—stock market
behavior is through network analysis, which can offer insights into the
relationships between stocks and the overall market structure. In this paper,
we seek to address the question: Can network analysis of the stock market,
specifically in observation of the lead-lag effect, provide valuable insights
for investors and market analysts? This inquiry is both interesting and
pertinent for several reasons.

Firstly, grasping the relationships between stocks and the overall market
structure can aid investors in making more informed, and potentially more
profitable, decisions regarding their investments. Additionally, network
analysis may offer new tools for monitoring the stock market and identifying
trends or potential risks. Through this investors will be able to observe the
bad returns of a leading stock and therefore draw conclusions about the
effects this will have on other closely tied securities.

To tackle this question, we will build upon two previous studies—of which both
observed a power-law distribution in stock returns. The first, performed by
researchers at HIT China, constructed a network for the US stock market based
on a model in which two stocks were said to be connected if their returns fell
in the respective returns of each other over some period to some threshold of
accuracy. The study found that through the use of this technique, they were
able to select a reasonably performing portfolio that beat the S&P 500 [1].
The second study, performed by a group of researchers at Iran’s Tehran
University, employed community detection techniques to construct a correlation
network. They discovered that the resulting communities were consistent with
market sectors classified using the Standard Industrial Classification code.
They also utilized network analysis to construct visualizations of the return
correlations among various public stocks, which offered an intuitive way to
examine the overall correlation structure of different public stocks and
identify key market segments [2].

While these prior works provide valuable insights into the network structure
of the stock market, there remains much to expand on this approach. In our
research, we concentrate on the lead-lag effect, which refers to the
phenomenon where the returns of one stock lead, over some period—referred to
as the "lag"—the returns of another stock. By analyzing the lead-lag effect
within the network of the stock market, we aim to offer insights into the
dynamics of market behavior and inform investment strategies.

##  2 Methodology

We realize that the lead-lag relationship is a "strategy-enhancer," but not a
strategy in itself. Thus, we combined it with a known, semi-reliable,
fundamental strategy: the Capital Asset Pricing Model (CAPM) as well as with
another network theory technique, out-degree. Simply put, CAPM is an effective
way to compute the expected return of a given security based on the market
risk premium and the systematic risk of the security and out-degree in a
measure of the number of edges leaving a given node. With out-degree,
specifically, we looked at the out-degree of the lagger—assuming that a lagger
with a lower out-degree will receive greater input from its few leaders. CAPM
is normalized around the market rate, so we curve our out-degree scores and
then weigh the two methods 70-30, CAPM and out-degree respectively. From this,
we choose the top six stock pairs that have the highest absolute sum, these
are the stocks that we will trade over the quarter. With the coupling of these
two strategies, we selected a portfolio and compared its returns to the S&P
500.

###  2.1 Data Collection and Preprocessing

We first constructed a network of the market using the stock return data of
the S&P 500. This involved building a script that web scraped for relevant
data. Using clever list manipulation and Python’s numba package, we were able
to tame the large computational time that comes with processing the data of
500 stocks for ten years. Correspondingly, there are roughly 250 trading days
per year so making sure that our data cubes matched imported S$\&$P 500 time
series data was critical to ensure a reliable comparison.

The data used in subsequent steps is sourced from a collections of files we
generated. One contains daily stock prices for all assets in the S$\&$P 500,
another contains the lead-lag preprocessed adjacency matrices in the form of a
tensor, and another contains quarterly CAPM data for all assets.

###  2.2 Identifying Lead-Lag Pairs

We then created our algorithm, selecting connection thresholds found in the
literature, first by constructing a 500x500 matrix of the closing prices of
the stocks. If, over consecutive time periods, the return of stock i stays
within a given interval of stock j then we consider stock j to lead stock i,
where the time period is called the lag. If a stock leads another stock we
placed a 1 in our adjacency matrix 111Which is a pseudo adjacency matrix, as
although it resembles an adjacency matrix it is not symmetric nor does it have
1’s completely along the diagonal.. Noticeably, stocks tend to lead
themselves, but because of volatility, certain stocks do not.

To identify the leader-lagger pairs, a 2D array is created from the summed
data, with the diagonal masked to ignore self-pairs, in other words, we ignore
stocks leading or lagging themselves. The top 20 largest values in the masked
array are then extracted, and the corresponding row and column indices are
identified. These indices are used to create a dictionary of leader-lagger
pairs, where the leader is the stock that leads in price movement, and the
lagger is the stock that follows. Using CAPM information of laggers, we update
this dictionary quarterly to select the top six best candidates222These pairs
have laggers with the lowest out-degree as mentioned in the Introduction.
according to our strategy and trade with them through the quarter.

###  2.3 Trading Strategy and Portfolio Simulation

The trading strategy used in this study is based on the following rules:

  * •

If a leader stock’s price increases by at least the buy threshold from the
previous day, buy the corresponding lagger stock.

  * •

If a leader stock’s price decreases below a trailing stop level, sell the
corresponding lagger stock.

The study assumes an initial investment of 500,000 USD, with the investment
amount evenly distributed among the lagger stocks. A portfolio dictionary is
initialized to track the number of shares held for each lagger stock. For each
trading day, the strategy’s rules are applied, and the portfolio is updated
quarterly with new CAPM data. Buying threshold and trailing stop percentage
are used in subsequent parameter exploration and used to maximize portfolio
return.

###  2.4 Performance Evaluation

The performance of the trading strategy is evaluated by calculating the
portfolio return over the period. The return is calculated as the difference
between the final value and the initial value of the portfolio, divided by the
initial value. The portfolio return is compared with the S$\&$P 500 index
return over the same period to provide context for the strategy’s performance
using the yfinance API. To visualize the performance, a line chart is created,
showing the cumulative returns for both the portfolio and the S$\&$P 500
index, see Fig. 2

###  2.5 Visualization of Leader-Lagger Relationships

To visualize the relationships between leader and lagger stocks, a directed
graph is created using Python’s networkx package, see Fig. 1. The nodes in the
graph represent stocks, while the directed edges represent the leader-lagger
relationships. Leader stocks are colored red, and lagger stocks are colored
blue. The directed graph is displayed using using Python’s matplotlib package.

##  3 Results

Our created portfolio beat the S$\&$P 500 by about ten percent over the course
of the last four fiscal quarters, ubiquitously known to be a bear market, see
Fig. 2(b).333We analyzed the lead-lag effect on data spanning from March 15,
2022 to March 15, 2023. Moreover, during a bull market, such as during the
second and third fiscal quarters of 2021, back-testing delivered a larger
positive portfolio return when compared to the S$\&$P 500 during the same
period, see Fig. 2(a). As well, in the major bull market of 2021 we beat the
S$\&$P 500 by about ten percent once again.

###  3.1 Network and Graph Model

Evident in Fig. 1(b), a network plot of our lead-lag pairs on the first day of
trading444It would be interesting to study the temporal evolution of these
plots, specifically the change in unique lead-lag pairs. during a bear market,
multiple leaders led the same lagger and the strongest signals often came from
one lagger. Interestingly enough, lead-lag pairs did not always form along
industry lines, in contradiction with previous research.555Likely because we
used a shorter lag period. Correspondingly, we observe a greater number of
lead-lag pairs during the first trading day of the bull market, see Fig. 1(a).
Since our analysis is centered around a directed graph, we are precluded from
studying the assortment of network centrality measures, namely, eigenvector
centrality.666In future work, the adjacency matrix can be made symmetric and
these measures can be deployed.

(a) Bull market lead-lag pairs.

(b) Bear market lead-lag pairs.

Figure 1: Network graphs of lead-lag pairs for the first day of trading. Blue
nodes represent laggers with directed edges connecting them to red nodes
representing leaders.

###  3.2 Portfolio Simulation

Over the course of the last four fiscal quarters and during the summer of
2021, our portfolio returns mirrored the ebb and flow of the S$\&$P 500’s
historical returns, see Fig. 2. Since bear markets are hard to predict it is
was reasonable for our portfolio simulation to not experience positive
returns777It is quite difficult without use of leverage to experience positive
gains in a bear market and our portfolio does not include an options trading
functionality., though we lose roughly five percent in value while the S$\&$P
500 approximately loses fifteen percent, see Fig. 2(b).888Hence, someone who
did not invest at all would have kept the value of their capital while our
portfolio and the S$\&$P 500’s value diminished, not accounting for inflation.
Moreover, during a bull market, we observe a positive return roughly twelve
percent greater than the S$\&$P 500 during the same period, see Fig. 2(a).

(a) Bull market of Q2 and Q3 of 2021.

(b) Bear market of the past four quarters.

Figure 2: Back-testing of portfolio returns compared to historical returns of
the S$\&$P 500 (∧GSPC).

###  3.3 Parameter Exploration

Naturally, we were prompted to conduct basic parameter exploration with the
hopes of maximizing our return given our risk preferences amongst other
variables.

####  3.3.1 Bear Market

We observed that when the buy threshold was set to a value greater than or
equal to 1.13, the model never engaged in any trades, thereby resulting in no
portfolio returns. This suggests that a buy threshold of 1.13 or higher is too
conservative for the model to capitalize on market opportunities.

Furthermore, we found that the highest portfolio value was attained when the
trailing stop percentage was set to a value less than or equal to 0.40. At
this threshold, all trades were sold on a daily basis, indicating a high
turnover strategy. This result is demonstrated in the contour plot presented
in Fig. 3(a).

In contrast, when the buy threshold was set to a value less than or equal to
1.07, the model traded more frequently, exposing the portfolio to increased
risk. The 2D cross-section of the contour plot in Fig. 3(c) shows the impact
of the buy threshold on trading frequency and risk exposure. As a result we
used a trailing stop percentage of 0.10 in our portfolio simulations and
toggled the buying threshold according to whether or not we were trading in a
bear or bull market, see Fig. 2.

(a) Contour plot for bear market parameters.

(b) 2D cross-section of bear contour plot, holding buy threshold constant.

(c) 2D cross-section of bear contour plot, holding stop loss threshold
constant.

Figure 3: Mesh grid plot of trailing stop percentage vs. buy threshold vs.
portfolio returns in a bear market.

####  3.3.2 Bull Market

The behavior observed in the bull market was, as expected, different to that
observed in a bear market.

We observed that when the buy threshold was set to a value greater than or
equal to 1.08999Although this seems contrary to intuition this is an artifact
of stocks in 2021 having lesser daily gains but but greater gains overall.
Moreover, 2022 had more visible return volatility but overall maintained this
behavior., the model refrained from engaging in any trades (Fig 4(c)),
consequently yielding no portfolio returns. This phenomenon may be attributed
to the fact that in a bull market, the model perceives a buy threshold of 1.08
or higher as overly conservative, thus missing out on potential market
opportunities. Additionally, we discovered that the trailing stop percentage
negatively impacted the portfolio when set above 0.80 (Fig 4(b)). This could
be due to the fact that a higher trailing stop in a bull market may lead to
premature exits from profitable positions, limiting the potential gains.

Interestingly, our study found that the optimal buy threshold was set at 1.02,
where the return increased from 1 to 1.02 and subsequently decreased up to
1.08( Fig 4(a)). This suggests that a buy threshold set within this range
allows the model to effectively capitalize on the upward market trend while
minimizing the risk of overexposure. The reason for this increase in return up
to 1.02 and subsequent decrease can be attributed to the balance between
capturing gains during the bull market and avoiding excessive risk exposure.
At a buy threshold of 1.02, the model is able to participate in the market’s
upward momentum, while still maintaining a relatively conservative approach to
minimize potential losses. As the buy threshold increases beyond 1.02, the
model becomes increasingly conservative, leading to missed opportunities and a
decline in returns.

(a) Contour plot for bull market parameters.

(b) 2D cross-section of bull contour plot, holding buy threshold constant.

(c) 2D cross-section of bull contour plot, holding stop loss threshold
constant.

Figure 4: Mesh grid plot of trailing stop percentage vs. buy threshold vs.
portfolio returns in a bull market.

In both bear and bull markets, our findings indicate that an optimal trading
strategy would involve setting the buy threshold and trailing stop percentage
within specific ranges to balance the trade-offs between portfolio returns,
trading frequency, and risk exposure. Further research is warranted to refine
these ranges and explore additional parameters that may influence the
performance of the model.

##  4 Analysis and Considerations

In our study on the lead-lag effect, there were several considerations that we
took into account while developing our model. Firstly, the lead-lag threshold
that we decided upon was not extensively tested, which could potentially
affect the robustness of our findings. There is also the issue of overfitting
when determining the buy threshold and stop loss parameters. Through parameter
estimation, we identified values that appeared to be profitable in both bull
and bear markets. However, it is crucial to exercise caution when generalizing
these findings to a larger sample, as overfitting can lead to overly
optimistic performance estimates that may not hold true in different market
conditions. Additionally, our chosen fundamental strategy, the Capital Asset
Pricing Model (CAPM), is based on certain unrealistic assumptions that may not
accurately reflect real-world trading conditions. Lastly, it is important to
note that due to the frequent trading nature of our model, commission costs
must be considered as they could significantly impact the overall performance
of the model. These considerations highlight the challenges and limitations
faced when attempting to develop an effective trading model that accurately
captures the complexities of the lead-lag effect.

##  5 Conclusions

Based on our analysis of the lead-lag effect in the NYSE, we draw several
significant conclusions. Firstly, our model, which combines network analysis
with fundamental strategies, can potentially outperform the market in
predicting stock performance. However, we acknowledge that our results are
based on a specific subset of years and may not be able to be generalized.

Secondly, our finding that lead-lag pairs often form outside of industry lines
contradicts previous literature and suggests that investors may benefit from
examining unconventional connections.

Thirdly, we demonstrate that our strategy can be customized to fit individual
investors’ risk preferences by adjusting the buy threshold. In bearish
markets, investors can raise the threshold to experience the market’s ups and
downs to a lesser degree, while in growth-favoring markets, investors can
lower the threshold to expose themselves to more risk and potential return.

Overall, our study highlights the usefulness of network analysis in predicting
stock performance and the importance of considering unconventional connections
in investment decisions. Furthermore, our approach allows for customization to
meet individual risk preferences, demonstrating its practical relevance to
investors.

##  6 Future Work

There are several potential areas for future work and unanswered questions
that could be addressed in further research. First, our study did not
thoroughly explore the impact of different time periods on the effectiveness
of the lead-lag effect. Investigating whether the relationships observed in
our study hold true over longer time periods or during different economic
cycles and market conditions would be an interesting avenue for future
research.

In addition, our study focused on the use of network analysis combined with
fundamental strategies in predicting stock performance. However, there are
several other strategies and factors that could be incorporated into the
analysis. For instance, investigating the effectiveness of incorporating
technical analysis, sentiment analysis, or macroeconomic indicators in
combination with network analysis could yield interesting results.

Moreover, it would be more realistic to incorporate weighted edges and a
distributed value for the lead-lag coefficient in our model, taking into
account the differences in market capitalization among stocks. Specifically,
assigning higher numbers in the matrix to stocks that were extremely lagged
over a given period, and assigning less weight to the edges of stocks with
higher market capitalization. As well, it would be helpful to investigate the
effectiveness of different lead-lag thresholds in predicting stock
performance.

## References

[1] Li, Y., Wang, T., Sun, B., & Liu. Detecting the lead–lag effect in stock
markets: Definition, patterns, and investment strategies. Financial
Innovation, 2022.

[2] A. Namaki, A.H. Shirazi, R. Raei, and G.R. Jafari. Network analysis of a
financial market based on genuine correlation and threshold method. Physica A:
Statistical Mechanics and its Applications, 2011.

[3] R. D. Smith. The Spread of the Credit Crisis: View from a Stock
Correlation Network. Journal of Korean Physical Society, 54:2460, June 2009.


11institutetext:  Institute of Mathematics, University of Kassel, Germany
11email<EMAIL_ADDRESS>
11email<EMAIL_ADDRESS>22institutetext: Faculty of
Mathematics, University of Vienna, Austria
22email<EMAIL_ADDRESS>
# A SageMath Package for
Elementary and Sign Vectors with
Applications to Chemical Reaction Networks

Marcus S. Aichmayr 11 0009-0008-8362-2644 Stefan Müller 22 0000-0002-3541-7856
Georg Regensburger 11 0000-0001-7735-3726

###### Abstract

We present our SageMath package elementary_vectors for computing elementary
and sign vectors of real subspaces. In this setting, elementary vectors are
support-minimal vectors that can be determined from maximal minors of a real
matrix representing a subspace. By applying the sign function, we obtain the
cocircuits of the corresponding oriented matroid, which in turn allow the
computation of all sign vectors of a real subspace.

As an application, we discuss sign vector conditions for existence and
uniqueness of complex-balanced equilibria of chemical reaction networks with
generalized mass-action kinetics. The conditions are formulated in terms of
sign vectors of two subspaces arising from the stoichiometric coefficients and
the kinetic orders of the reactions. We discuss how these conditions can be
checked algorithmically, and we demonstrate the functionality of our package
sign_vector_conditions in several examples.

###### Keywords:

elementary vectors sign vectors oriented matroids generalized mass-action
systems deficiency zero theorem robustness.

††footnotetext:  The SageMath packages are available at:
https://github.com/MarcusAichmayr/
The examples of this document are available at:
https://marcusaichmayr.github.io/sign_vector_conditions/

##  1 Elementary vectors of a subspace

For real subspaces, elementary vectors are nonzero vectors with minimal
support, as introduced in [13]. Since we also deal with parameters, we
consider elementary vectors for subspaces $Q^{n}$ of the quotient field $Q$ of
an integral domain $R$. The _support_ of a vector $x\in Q^{n}$ is the index
set of nonzero components,

$\operatorname{supp}x=\\{i\mid x_{i}\neq 0\\}.$

###### Definition 1

For a subspace $\mathcal{V}$ of $Q^{n}$, we call a nonzero vector
$v\in\mathcal{V}$ _elementary_ if, for all nonzero vectors $x\in\mathcal{V}$,

$\operatorname{supp}x\subseteq\operatorname{supp}v\quad\text{implies}\quad\operatorname{supp}x=\operatorname{supp}v.$

Elementary vectors with the same support are easily seen to be multiples.
Therefore, a subspace contains only finitely many elementary vectors up to
multiples. Further, a subspace is always generated by a finite set of
elementary vectors.

We represent a subspace as the kernel of a matrix $M\in R^{d\times n}$ with
rank $d$ over $Q$. If the kernel of a matrix is $1$-dimensional, it has
exactly one elementary vector. We apply this fact by considering submatrices
consisting of $d+1$ columns. To find kernel vectors, we solve systems of
linear equations. In particular, we compute maximal minors and apply Cramer’s
Rule. By inserting additional zeros, we extend the kernel vectors to kernel
vectors of the initial matrix. The resulting vectors are elementary if they
are nonzero.

For $I\subseteq[n]=\\{1,\ldots,n\\}$, we denote by $M_{I}$ the submatrix of
$M$ that consists only of the columns corresponding to the indices $I$. We
obtain a formula for computing elementary vectors in $R^{n}$.

######  Proposition 1 (cf. equation (2.1) in [6])

For a matrix $M\in R^{d\times n}$ with rank $d$ and $I\subseteq[n]$ with
$\lvert I\rvert=d+1$, define the vector $v\in R^{n}$ by

$v_{i}=\begin{cases}(-1)^{\lvert\\{k\in I\colon k<i\\}\rvert}\det
M_{I\setminus\\{i\\}},&\text{if $i\in I$},\\\ 0,&\text{otherwise}.\end{cases}$

The vector $v\in\ker M$ is elementary if $\operatorname{rank}M_{I}=d$.

For computing all elementary vectors in $\ker M$, we need to compute
$\binom{n}{d}$ determinants of all $d\times d$ submatrices of $M$ (maximal
minors). In contrast, for Gaussian Elimination, we would need to compute the
kernel of a $(d+1)\times d$ matrix for each of the $\binom{n}{d+1}$ index sets
$I$. In our implementation, it turns out that the approach using maximal
minors is more efficient for computing elementary vectors for medium-sized
matrices.

If all maximal minors are nonzero, there is exactly one elementary vector for
each index set $I$. This gives an upper bound of $\binom{n}{d+1}$ elementary
vectors with pairwise distinct support of an $(n-d)$-dimensional subspace. If
the rank of $M_{I}$ is not maximal, we obtain the zero vector. Proposition 1
suggests an algorithm for computing elementary vectors, which we have
implemented in our SageMath package elementary_vectors [1]. We demonstrate it
by an example.

⬇

sage: from elementary_vectors import *

e: M = matrix([[1, 1, 2, 0], [0, 0, 1, 2]])

e: M.minors(2)

1, 2, 1, 2, 4]

e: elementary_vectors(M)

, -1, 0, 0), (4, 0, -2, 1), (0, 4, -2, 1)]

Note that the first maximal minor is zero. This is the reason why we obtain
only $3$ and not $\binom{4}{3}=4$ elementary vectors.

###  1.0.1 Solvability of linear inequality systems.

A fundamental theorem for deciding the solvability of linear inequality
systems is given in [13] and [14]. Note that every linear inequality system
can be written as an intersection of a subspace and a Cartesian product of
intervals. By iterating over elementary vectors, we check whether such an
intersection is empty. For more details, we refer to our manuscript [3] (see
also [7]).

######  Theorem 1.1 (“Minty’s Lemma”, Theorem 22.6 in [14])

For a subspace $\mathcal{V}$ of $\mathbb{R}^{n}$ and nonempty intervals
$\mathcal{I}=I_{1}\times\cdots\times I_{n}$,

  * either there exists a vector $x\in\mathcal{V}\cap\mathcal{I}$,

  * or there exists an elementary vector $v\in\mathcal{V}^{\perp}$ with $v^{\top}z>0$ for all $z\in\mathcal{I}$.

Based on Theorem 1.1, the function exists_vector decides the solvability of
such systems.

⬇

sage: from vectors_in_intervals import *

e: M = matrix([[1, 0, 1, 0], [0, 1, 1, 1]])

e: I = intervals_from_bounds([2,5,0,-oo],[5,oo,8,5],

.: [True,True,False,False],[False,False,False,True])

e: I

, 5), [5, +oo), (0, 8), (-oo, 5]]

e: exists_vector(M, I)

e

###  1.0.2 Sign vectors.

We call elements in $\\{-,0,+\\}^{n}$ _sign vectors_. For
$x\in\mathbb{R}^{n}$, we define the sign vector
$\operatorname{sign}(x)\in\\{-,0,+\\}^{n}$ by applying the sign function
componentwise. For a set $S\subseteq\mathbb{R}^{n}$, we obtain the set of sign
vectors

$\operatorname{sign}(S)=\\{\operatorname{sign}(x)\mid x\in
S\\}\subseteq\\{-,0,+\\}^{n}$

by applying sign to each vector of $S$. As for real vectors, we introduce the
_support_ $\operatorname{supp}\sigma=\\{i\mid\sigma_{i}\neq 0\\}$ of a sign
vector $\sigma$.

Further, we obtain a partial order on $\\{-,0,+\\}^{n}$ by defining $0<-,+$.
As in [9], the _(lower) closure_ of a set of sign vectors
$T\subseteq\\{-,0,+\\}^{n}$ is the set

$\displaystyle\overline{T}=\\{\sigma\in\\{-,0,+\\}^{n}\mid\text{$\sigma\leq\tau$
for some $\tau\in T$}\\}.$

A (realizable) _oriented matroid_ is the set of sign vectors that correspond
to a real subspace. We call the elements of an oriented matroid _covectors_.
To obtain them, we apply the sign function to the elements in a subspace.

The sign vectors corresponding to the elementary vectors are called
_cocircuits_. They generate all elements of an oriented matroid, just like the
elementary vectors generate all elements of a subspace. Since the cocircuits
are determined by the signs of the maximal minors, called _chirotopes_ , we
can express many sign vector conditions in terms of chirotopes. For further
details on oriented matroids, we refer to [4, Chapter 7], [15, Chapters 2 and
6], [12] and the encyclopedic study [5]. Our package also offers several
functions for oriented matroids. As an example, we show the computation of
cocircuits.

⬇

sage: from sign_vectors.oriented_matroids import *

e: M = matrix([[1, 1, 2, 0], [0, 0, 1, 2]])

e: cocircuits_from_matrix(M)

+00), (-0+-), (+0-+), (+-00), (0-+-), (0+-+)}

e: covectors_from_matrix(M)

000), (-+00), (--+-), (+-00), (+--+), (-0+-), (+0-+), (0-+-), (0+-+), (++-+),
(-+-+), (+-+-), (-++-)}

##  2 Applications to chemical reaction networks

For chemical reaction networks (CRNs) with generalized mass-action kinetics,
we recall basic notions from [11]. See also [9] and [10].

A generalized mass-action system $(G_{k},y,\widetilde{y})$ is given by a
simple directed graph $G=(V,E)$, positive edge labels
$k\in\mathbb{R}^{E}_{>}$, and two maps ${y\colon V\to\mathbb{R}^{n}}$ and
${\widetilde{y}\colon V_{s}\to\mathbb{R}^{n}}$, where $V_{s}=\\{i\mid i\to
i^{\prime}\in E\\}\subseteq V$ denotes the set of source vertices. Every
vertex $i\in V$ is labeled with a (stoichiometric) complex
$y(i)\in\mathbb{R}^{n}$, and every source vertex $i\in V_{s}$ is labeled with
a kinetic-order complex $\widetilde{y}(i)\in\mathbb{R}^{n}$. Further, every
edge $(i\to i^{\prime})\in E$ is labeled with a rate constant $k_{i\to
i^{\prime}}>0$ and represents the chemical reaction ${y(i)\to y(i^{\prime})}$.
If every component of $G$ is strongly connected, $G$ and
$(G_{k},y,\widetilde{y})$ are called weakly reversible. The associated ODE
system for the positive concentrations $x\in\mathbb{R}^{n}_{>}$ (of $n$
chemical species) is given by

$\frac{\mathrm{d}x}{\mathrm{d}t}=\sum_{(i\to i^{\prime})\in E}k_{i\to
i^{\prime}}\,x^{\widetilde{y}(i)}\big{(}y(i^{\prime})-y(i)\big{)}.$ (1)

The sum ranges over all reactions, and every summand is a product of the
reaction rate $k_{i\to i^{\prime}}\,x^{\widetilde{y}(i)}$, involving a
monomial $x^{\widetilde{y}}=\prod_{j=1}^{n}(x_{j})^{\widetilde{y}_{j}}$
determined by the kinetic-order complex of the educt, and the reaction vector
$y(i^{\prime})-y(i)$ given by the stoichiometric complexes of product and
educt.

Let $I_{E},I_{E,s}\in\mathbb{R}^{V\times E}$ be the incidence and source
matrices of the digraph $G$, respectively, and

$A_{k}=I_{E}\operatorname{diag}(k)(I_{E,s})^{\top}\in\mathbb{R}^{V\times V}$

be the Laplacian matrix of the labeled digraph $G_{k}$. (This definition is
used in dynamical systems. In other fields, the Laplacian matrix is defined as
$A_{k}^{\top}$, $-A_{k}$, or $-A_{k}^{\top}$.) Now, the right-hand-side of (1)
can be decomposed into stoichiometric, graphical, and kinetic-order
contributions,

$\frac{\mathrm{d}x}{\mathrm{d}t}=YI_{E}\operatorname{diag}(k)(I_{E,s})^{\top}\,x^{\widetilde{Y}}=YA_{k}\,x^{\widetilde{Y}},$
(2)

where $Y,\,\widetilde{Y}\in\mathbb{R}^{n\times V}$ are the matrices of
stoichiometric and kinetic-order complexes, and
$x^{\widetilde{Y}}\in\mathbb{R}^{V}_{>}$ denotes the vector of monomials, that
is, $(x^{\widetilde{Y}})_{i}=x^{\widetilde{y}(i)}$. Clearly, the change over
time lies in the stoichiometric subspace $S=\operatorname{im}(YI_{E}),$ that
is, $\frac{\mathrm{d}x}{\mathrm{d}t}\in S$. Equivalently, trajectories are
confined to cosets of $S$, that is, $x(t)\in x(0)+S$. For positive
$x^{\prime}\in\mathbb{R}^{n}_{>}$, the set
$(x^{\prime}+S)\cap\mathbb{R}^{n}_{>}$ is called a stoichiometric class.

The (stoichiometric) deficiency is given by

$\delta=\dim(\ker Y\cap\operatorname{im}I_{E})=|V|-\ell-\dim(S),$

where $|V|$ is the number of vertices, and $\ell$ is the number of connected
components of the digraph. Analogously, we introduce the kinetic-order
subspace $\widetilde{S}=\operatorname{im}(\widetilde{Y}I_{E})$ and the
kinetic(-order) deficiency
$\widetilde{\delta}=\dim(\ker\widetilde{Y}\cap\operatorname{im}I_{E})=|V|-\ell-\dim(\widetilde{S}).$

A steady state $x\in\mathbb{R}^{n}_{>}$ of (2) that fulfills
$A_{k}\,x^{\widetilde{Y}}=0$ is called a positive complex-balanced equilibrium
(CBE). On the one hand, if $\delta=0$, then every equilibrium is complex-
balanced. On the other hand, if there exists a CBE, then the underlying graph
is weakly reversible, see the comment at the end of Section 4 in [11], and cf.
[10, Proposition 2.18].

######  Theorem 2.1 (robust $\delta=\widetilde{\delta}=0$ theorem, Theorem 46
in [9])

For a generalized mass-action system, there exists a unique positive CBE in
every stoichiometric class, for all rate constants and for all small
perturbations of the kinetic orders, if and only if
$\delta=\widetilde{\delta}=0$, the network is weakly reversible, and
$\operatorname{sign}(S)\subseteq\overline{\operatorname{sign}(\widetilde{S})}$.

We consider a CRN given by a graph with 5 vertices and 6 edges in 2 connected
components, and we label the vertices with stoichiometric and kinetic-order
complexes, respectively. (The two resulting labeled graphs are shown below.)
The kinetic-order complexes involve parameters $a,b,c\in\mathbb{R}$.

$A+B$$C$$D$$A$$E$$aA+bB$$C$$cA+D$$A$$E$

The resulting stoichiometric and kinetic-order subspaces are given by

$S=\operatorname{im}\begin{array}[]{cc}{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>B\\\ C\\\ D\\\ E\end{array}}&\begin{pmatrix}-1&0&-1\\\ -1&0&0\\\ 1&-1&0\\\
0&1&0\\\
0&0&1\end{pmatrix}\end{array}\quad\text{and}\quad\widetilde{S}=\operatorname{im}\begin{pmatrix}-a&c&-1\\\
-b&0&0\\\ 1&-1&0\\\ 0&1&0\\\ 0&0&1\end{pmatrix}.$

Clearly, the network is weakly reversible, and it is easy to verify that
$\delta=\widetilde{\delta}=0$. To study existence and uniqueness of complex-
balanced equilibria, we compute the sign vectors of $S$ and $\widetilde{S}$.

⬇

sage: from sign_vectors.oriented_matroids import *

e: S = matrix([[-1,-1,1,0,0],[0,0,-1,1,0],[-1,0,0,0,1]])

e: covectors_from_matrix(S, kernel=False)

0000), (00+-0), (0-0+-), (+000-), (+-++-), (0--+-), (+0+--), (+++--), (0-+--),
(---+-), (--+0-), (--+--), (++0-0), ..., (++-++), (--+0+), (---++), (-+0-+)}

Since $\widetilde{S}$ depends on the parameters $a,b$ and $c$, the sign
vectors also depend on these parameters and hence, we cannot compute them
directly. Thus, we consider specific values for the parameters and determine
the corresponding sign vectors.

⬇

sage: var(’a, b, c’);

e: St = matrix([[-a,-b,1,0,0],[c,0,-1,1,0],[-1,0,0,0,1]])

e: covectors_from_matrix(St(a=2, b=1, c=1), kernel=False)

0000), (+-++-), (0--+-), (+0+--), (+++--), (0-0+-), (0-+--), (---+-), (00-++),
(+000-), (+0-+0), (--+--), (0++-0), ..., (++-+0), (-0+-0), (-+0-+), (-0+--)}

For $a=2$, $b=1$ and $c=1$,
$\operatorname{sign}(S)\subseteq\overline{\operatorname{sign}(\widetilde{S})}$.
Consequently, the conditions of Theorem 2.1 are satisfied in this specific
case. Obviously, we cannot cover all possible cases for the parameters that
way. However, by expressing this condition in terms of maximal minors of the
kernel matrices, we can compute with the parameters directly.

######  Proposition 2 (Proposition 32 in [9])

For subspaces $S$, $\widetilde{S}\subseteq\mathbb{R}^{n}$ of dimension $n-d$
and matrices $W,\,\widetilde{W}\in\mathbb{R}^{d\times n}$ with $S=\ker W$,
$\widetilde{S}=\ker\widetilde{W}$, and rank $d$, the following are equivalent:

  1. 1.

$\operatorname{sign}(S)\subseteq\overline{\operatorname{sign}(\widetilde{S})}$.

  2. 2.

$\det W_{I}\neq 0$ implies $\det W_{I}\det\widetilde{W}_{I}>0$ for all subsets
$I\subseteq[n]$ with $\lvert I\rvert=d$ (or “$<0$” for all $I$).

In the example, $S=\ker W$ and $\widetilde{S}=\ker\widetilde{W}$ with

$\displaystyle W=\begin{pmatrix}1&0&1&1&1\\\
0&1&1&1&0\end{pmatrix}\qquad\text{and}\qquad\widetilde{W}=\begin{pmatrix}1&0&a&a-c&1\\\
0&1&b&b&0\end{pmatrix}.$

By applying Proposition 2, we obtain several conditions on $a$, $b$ and $c$.
We use package [2] to compute this (and several other) sign vector
condition(s).

⬇

sage: W = matrix([[1, 0, 1, 1, 1], [0, 1, 1, 1, 0]])

e: var(’a, b, c’);

e: Wt = matrix([[1, 0, a, a - c, 1], [0, 1, b, b, 0]])

e: from sign_vector_conditions import *

e: condition_closure_minors(W, Wt)

- c > 0, b > 0, a > 0}]

Hence, the network has a unique positive CBE if and only if $a,b>0$ and $a>c$.

###  2.0.1 Uniqueness of CBE.

Also the uniqueness of CBE (in every stoichiometric class and for all rate
constants) — and hence its converse: multiple CBE – can be characterized in
terms of a sign vector condition. For further details on sign vector
conditions for injectivity in the context of CRNs and further references, we
refer to [8].

######  Proposition 3 (cf. Proposition 3.1 in [10])

For a generalized mass-action system, there exists at most one positive CBE in
every stoichiometric class, for all rate constants, if and only if

$\operatorname{sign}(S)\cap\operatorname{sign}(\widetilde{S}^{\perp})=\\{0\\}.$
(3)

Because of the parameters, we cannot directly compute the sign vectors in the
kernel of $\widetilde{S}$. Again, we use maximal minors to express (3).

######  Corollary 1 (cf. Corollary 4 in [9])

For subspaces $S$, $\widetilde{S}\subseteq\mathbb{R}^{n}$ of dimension $n-d$
and matrices $W,\,\widetilde{W}\in\mathbb{R}^{d\times n}$ with $S=\ker W$,
$\widetilde{S}=\ker\widetilde{W}$, and rank $d$, the following are equivalent:

  1. 1.

$\operatorname{sign}(S)\cap\operatorname{sign}(\widetilde{S}^{\perp})=\\{0\\}$.

  2. 2.

Either $\det W_{I}\det\widetilde{W}_{I}\geq 0$ for all $I\subseteq[n]$ with
$\lvert I\rvert=d$,
or $\det W_{I}\det\widetilde{W}_{I}\leq 0$ for all $I\subseteq[n]$ with
$\lvert I\rvert=d$.

Comparing the maximal minors yields:

⬇

sage: condition_uniqueness_minors(W, Wt)

- c >= 0, a >= 0, b >= 0}]

Hence, positive CBE are unique if and only if $a,b\geq 0$ and $a\geq c$.

###  2.0.2 Unique existence of CBE.

We discuss a novel algorithm for checking a certain degeneracy condition for
subspaces that is part of a characterization of the unique existence of CBE.
The result also involves the set of nonnegative sign vectors
$T_{\oplus}=T\cap\\{0,+\\}^{n}$ of a set of sign vectors
$T\subseteq\\{-,0,+\\}^{n}$.

######  Theorem 2.2 ($\delta=\widetilde{\delta}=0$ theorem, Theorem 45 in [9])

For a generalized mass-action system, there exists a unique positive CBE in
every stoichiometric class, for all rate constants, if and only if
$\delta=\widetilde{\delta}=0$, the network is weakly reversible, and

  1. 1.

$\operatorname{sign}(S)\cap\operatorname{sign}(\widetilde{S}^{\perp})=\\{0\\}$;

  2. 2.

for all nonzero
$\widetilde{\tau}\in\operatorname{sign}(\widetilde{S}^{\perp})_{\oplus}$,
there is a nonzero $\tau\in\operatorname{sign}(S^{\perp})_{\oplus}$ with
$\tau\leq\widetilde{\tau}$; and

  3. 3.

$(S,\widetilde{S})$ is nondegenerate.

We discussed the first condition above. To check the second condition, we use
nonnegative cocircuits. For the third condition, we reformulate Definition 13
in [9], regarding the (non-)degeneracy of two subspaces.

###### Definition 2

A pair $(S,\widetilde{S})$ of subspaces of $\mathbb{R}^{n}$ is called
_degenerate_ if there exists $z\in\widetilde{S}^{\perp}$ with a positive
component such that

  1. 1.

for all $I_{\lambda}=\\{i\in[n]\mid z_{i}=\lambda\\}$ for some $\lambda>0$,
there exists $\pi\in\operatorname{sign}(S)_{\oplus}$ such that $\pi_{i}=+$ iff
$i\in I_{\lambda}$; and

  2. 2.

for all nonzero $\tau\in\operatorname{sign}(S^{\perp})_{\oplus}$, we have
$\operatorname{supp}\tau\not\subseteq\operatorname{supp}z$.

That is, the pair $(S,\widetilde{S})$ is degenerate if there exists
$z\in\widetilde{S}^{\perp}$ such that its equal positive components are
covered by nonnegative sign vectors of $S$, and there is no nonzero,
nonnegative sign vector in $S^{\perp}$ such that its support is contained in
$\operatorname{supp}z$.

Since we simply cannot iterate over the subspace $\widetilde{S}^{\perp}$, we
have to find a different approach. Instead, we consider sets of nonnegative
covectors in $\operatorname{sign}(S)$. Then, we use Theorem 1.1 (Minty’s
Lemma) to decide whether a vector $z$ exists that has positive equal
components on the support of each of these covectors.

We reformulate Definition 2 by demanding the existence of a set of nonnegative
covectors that cover all equal positive components of a vector $z$.

###### Definition 3

A pair $(S,\widetilde{S})$ of subspaces in $\mathbb{R}^{n}$ is _degenerate_ if
there exists a set of nonzero covectors
$T\subseteq\operatorname{sign}(S)_{\oplus}$ with disjoint support and a vector
$z\in\widetilde{S}^{\perp}$ such that

  1. 1.

    1. (a)

for all $\pi\in T$, $\lambda_{\pi}=z_{i}=z_{j}>0$ for all
$i,j\in\operatorname{supp}\pi$,

    2. (b)

for all $i\notin\bigcup_{\pi\in T}\operatorname{supp}\pi$, $z_{i}\leq 0$, and

    3. (c)

the $\lambda_{\pi}$’s are pairwise distinct;

and

  2. 2.

for all nonzero $\tau\in\operatorname{sign}(S^{\perp})_{\oplus}$, we have
$\operatorname{supp}\tau\not\subseteq\operatorname{supp}z$.

If two $\lambda_{\pi}$’s are equal, we compose the corresponding covectors to
cover equal components of $z$. Therefore, Condition (1c) is redundant. For the
same reason, the supports of the covectors do not need to be disjoint. Since
nonnegative covectors can be represented as a composition of nonnegative
cocircuits, it suffices to consider cocircuits instead of covectors. Further,
note that we can check Condition 2 using cocircuits. Following these
observations, we obtain another reformulation of Definition 2.

###### Definition 4

A pair $(S,\widetilde{S})$ of subspaces in $\mathbb{R}^{n}$ is _degenerate_ if
there exists a set of nonnegative cocircuits
$C\subseteq\operatorname{sign}(S)_{\oplus}$ and a vector
$z\in\widetilde{S}^{\perp}$ such that

  1. 1.

    1. (a)

for all $\pi\in C$, $z_{i}=z_{j}>0$ for all $i,j\in\operatorname{supp}\pi$,
and

    2. (b)

for all $i\notin\bigcup_{\pi\in C}\operatorname{supp}\pi$, $z_{i}\leq 0$;

and

  2. 2.

for all cocircuit $\tau\in\operatorname{sign}(S^{\perp})_{\oplus}$, we have
$\operatorname{supp}\tau\not\subseteq\operatorname{supp}z$.

To check degeneracy algorithmically, we iterate over sets of nonnegative
cocircuits. In particular, Algorithm 1 below recursively constructs such sets
and determines whether a corresponding $z$ exists. One could store the
cocircuits and construct a subspace of $\widetilde{S}^{\perp}$ using the
conditions on equal components. Here, we modify the subspace such that its
elements are equal on the components corresponding to the cocircuits in each
step (line 1) and keep track of the positive entries using an index set. If
this subspace contains a vector that is positive on exactly this index set
(line 1), condition 1 holds. If this vector also satisfies condition 2 (line
1), it certifies degeneracy. For efficiency, we use cocircuits to check this
condition. Note that we apply exists_vector from Section 1.0.1 to efficiently
check for existence in line 1 and 1.

1 Function  _degenerate(_$S,\widetilde{S}$_)_:

2 $C:=$ set of nonnegative cocircuits of $\operatorname{sign}(S)$

3 $\mathcal{V}:=\widetilde{S}^{\perp}$

4 global $\texttt{{is_degenerate}}:=\texttt{{False}}$

5 recursive(_$C,\mathcal{V},\emptyset$_)

6 return is_degenerate

7

8Function  _recursive(_$C,\mathcal{V},I$_)_:

9 while  _$C\neq\emptyset$_ do

10 choose any $\pi\in C$

11 $C:=C\setminus\\{\pi\\}$

12 $\overline{\mathcal{V}}:=$ subspace of $\mathcal{V}$ where vectors are
equal on $\operatorname{supp}\pi$

13 $\overline{I}:=I\cup\operatorname{supp}\pi$

14 if  _$z\in\overline{\mathcal{V}}$ exists with $z_{i}>0$ iff
$i\in\overline{I}$_ then

15 for  _$\sigma\in\operatorname{sign}(\overline{\mathcal{V}})$ with
$\sigma_{i}=+$ iff $i\in\overline{I}$_ do

16 if  _ $\operatorname{supp}\tau\not\subseteq\operatorname{supp}\sigma$ for
all cocircuits $\tau\in\operatorname{sign}(S^{\perp})_{\oplus}$ _ then

17 $\texttt{{is_degenerate}}:=\texttt{{True}}$

18 return

19

20

21

22 else if  _$z\in\overline{\mathcal{V}}$ exists with $z_{i}>0$ if
$i\in\overline{I}$_ then

23 recursive(_$C,\overline{\mathcal{V}},\overline{I}$_)

24

25 if  _is_degenerate_ then return

26

27 return

28

Algorithm 1 subspaces degenerate

We consider Example 20 from [9]. Here, we have matrices

$\displaystyle W=\begin{pmatrix}0&0&1&1&-1&0\\\ 1&-1&0&0&0&-1\\\
0&0&1&-1&0&0\end{pmatrix}\quad\text{and}\quad\widetilde{W}=\begin{pmatrix}1&1&0&0&-1&a\\\
1&-1&0&0&0&0\\\ 0&0&1&-1&0&0\end{pmatrix}.$

The existence of a unique positive CBE depends on $a>0$.

⬇

sage: var(’a’); assume(a > 0);

e: W=matrix(3,6,[0,0,1,1,-1,0,1,-1,0,0,0,-1,0,0,1,-1,0,0])

e: Wt=matrix(3,6,[1,1,0,0,-1,a,1,-1,0,0,0,0,0,0,1,-1,0,0])

The first two conditions of Theorem 2.2 are independent of $a$.

⬇

sage: condition_uniqueness_sign_vectors(W, Wt)

e

e: condition_faces(W, Wt)

e

Condition 3 holds iff $a\in(0,1)\cup(1,2)$ as we demonstrate for specific
values.

⬇

sage: condition_nondegenerate(W, Wt(a=1/2))

e

e: condition_nondegenerate(W, Wt(a=2))

se

###  2.0.3 Acknowledgements

This research was funded in part by the Austrian Science Fund (FWF), grant
10.55776/P33218 to SM.

###  2.0.4

The authors have no competing interests to declare.

## References

  * [1] Aichmayr, M.S.: elementary_vectors, https://github.com/MarcusAichmayr/elementary_vectors
  * [2] Aichmayr, M.S.: sign_vector_conditions, https://github.com/MarcusAichmayr/sign_vector_conditions
  * [3] Aichmayr, M.S., Regensburger, G.: How to certify solvability of linear inequality systems with elementary vectors (2024), in preparation
  * [4] Bachem, A., Kern, W.: Linear programming duality: an introduction to oriented matroids. Berlin: Springer-Verlag (1992)
  * [5] Björner, A., Las Vergnas, M., Sturmfels, B., White, N., Ziegler, G.: Oriented matroids. Cambridge: Cambridge University Press, 2nd ed. edn. (1999). https://doi.org/10.1017/CBO9780511586507
  * [6] Brualdi, R.A., Friedland, S., Pothen, A.: The sparse basis problem and multilinear algebra. SIAM Journal on Matrix Analysis and Applications 16, 1–20 (1995)
  * [7] Minty, G.J.: A ’from scratch’ proof of a theorem of Rockafellar and Fulkerson. Mathematical Programming 7, 368–375 (1974)
  * [8] Müller, S., Feliu, E., Regensburger, G., Conradi, C., Shiu, A., Dickenstein, A.: Sign conditions for injectivity of generalized polynomial maps with applications to chemical reaction networks and real algebraic geometry. Foundations of Computational Mathematics 16, 69–97 (2016). https://doi.org/10.1007/s10208-014-9239-3
  * [9] Müller, S., Hofbauer, J., Regensburger, G.: On the bijectivity of families of exponential/generalized polynomial maps. SIAM Journal on Applied Algebra and Geometry 3, 412–438 (2019). https://doi.org/10.1137/18M1178153
  * [10] Müller, S., Regensburger, G.: Generalized mass action systems: complex balancing equilibria and sign vectors of the stoichiometric and kinetic-order subspaces. SIAM Journal on Applied Mathematics 72, 1926–1947 (2012). https://doi.org/10.1137/110847056
  * [11] Müller, S., Regensburger, G.: Generalized mass-action systems and positive solutions of polynomial equations with real and symbolic exponents (invited talk). In: Computer algebra in scientific computing. 16th international workshop, CASC 2014. Proceedings, pp. 302–323. Berlin: Springer (2014). https://doi.org/10.1007/978-3-319-10515-4_22
  * [12] Richter-Gebert, J., Ziegler, G.M.: Oriented matroids. In: Handbook of discrete and computational geometry, pp. 111–132. Boca Raton, FL: CRC Press (1997)
  * [13] Rockafellar, R.T.: The elementary vectors of a subspace of $\mathbb{R}^{n}$. Combin. Math. Appl., Proc. Conf. Univ. North Carolina 1967, 104–127 (1969)
  * [14] Rockafellar, R.T.: Convex analysis. Princeton University Press, Princeton, NJ (1970)
  * [15] Ziegler, G.M.: Lectures on polytopes. Springer (1995). https://doi.org/10.1007/978-1-4613-8431-1