# Estimating Stable Fixed Points and Langevin Potentials for Financial
Dynamics

Tobias Wand  Corresponding author’s e-mail address:
<EMAIL_ADDRESS>Institut für Theoretische Physik, Westfälische
Wilhelms-Universität Münster  Center for Nonlinear Science Münster,
Westfälische Wilhelms-Universität Münster  Timo Wiedemann  Jan Harren  Finance
Center Münster, Westfälische Wilhelms-Universität Münster  Oliver Kamps
Center for Nonlinear Science Münster, Westfälische Wilhelms-Universität
Münster

###### Abstract

The Geometric Brownian Motion (GBM) is a standard model in quantitative
finance, but the potential function of its stochastic differential equation
(SDE) cannot include stable nonzero prices. This article generalises the GBM
to an SDE with polynomial drift of order $q$ and shows via model selection
that $q=2$ is most frequently the optimal model to describe the data.
Moreover, Markov chain Monte Carlo ensembles of the accompanying potential
functions show a clear and pronounced potential well, indicating the existence
of a stable price.

Langevin Equation; Stochastic Differential Equation; Finance, Econophysics;
Data-Driven Inference

††preprint: APS/123-QED

##  I Introduction

Empirical research on financial data with methods from physics has lead to a
better understanding of correlation matrices in financial markets [1, 2, 3],
uncovered scaling behaviour of empirical distributions [4] and validated
microscopic models of the behaviour of individual traders [5, 6].
Collectively, these efforts became known as the discipline of Econophysics and
have sparked numerous further research projects [7, 8].

Ever since Bachelier’s seminal work introduced the Brownian Motion to describe
stock prices and the expansion to the Geometric Brownian Motion (GBM) in the
Black-Scholes-Merton model [9, 10, 11], differential equations have been an
important tool to analyse financial data. Econophysicists have used
(stochastic) differential equations to analyse currency exchange data [12],
interest rates [13] and stock market crashes [14] or derived (partial)
differential equations from microscopic trading models [15, 16]. A recent
empirical study tried to model price time series with a harmonic oscillator
ODE to reconcile the randomness of financial markets with the idea of a fair
price [17]. The Geometric Brownian Motion, still widely used as a standard
model for price time series, presents the researchers with a subtle difficulty
with regards to its interpretation: Its deterministic part implies either an
unlimited exponential growth or an exponential decline to a price of 0 as
pointed out in [18, 19, 20]. While traditional finance models have tried to
improve the Geometric Brownian Motion by changing its stochastic component,
the deterministic part has largely been left unchanged (cf. the discussion in
section 1 of [20]). While [20] used a constrained model with regularisation
via strong prior information to fit parameters to their model, the goal of the
present article is to estimate model parameters and to select the best model
without any of these restrictions, i.e. letting the data speak for itself.

The estimation of Langevin equations from data via the Kramers-Moyal
coefficients [21, 22] sparked a family of methods to estimate nonparametric
drift and diffusion coefficients to model the observed system as a stochastic
differential equation which have also been applied to financial data [23]. A
particularly interesting expansion of this method is given by the maximum-
likelihood-framework (ML) in [24]: for each time step $t_{i}$ and observed
data $x_{i}$, the transition likelihood $L_{i}=p(x_{i+1}|x_{i})$ from $x_{i}$
to $x_{i+1}$ is calculated and the joint likelihood $L=\sum_{i}L_{i}$ is
maximised by the estimation algorithm. This approach takes the inherent
stochasticity of stochastic differential equations into account and can be
performed with a parametric model to recover algebraic equations increase its
interpretability. Similarly, the SINDy algorithm recovers a sparse functional
form of the underlying algebraic equations by fitting the data to a candidate
function library, but struggles with noisy and stochastic data [25, 26].

This article uses a combination of the ML framework of [24] with the candidate
function library in [25] for a robust method to estimate stochastic
differential equations from data similar to [27]. As an extension, the
presented method can be used to estimate data from time series with non-
constant time increments $dt_{i}\neq dt_{j}$. We use stock market prices at
daily and 30-minute intervals as described in section II to estimate their
stochastic differential equations. In particular, we estimate the potential in
which the dynamics take place to evaluate the stability of the dynamical
process as explained in section III. The results for the different polynomial
orders of the model and their implication for the stability are shown in
section IV and discussed in section V.

##  II Data

We analyse stock market data from the companies listed in table 1 to cover a
range of different business sectors. Our analysis covers two distinct market
conditions: (i) a calm period from early 2019 through early 2020 which was
characterised by low overall volatility and (ii) the Covid Selloff beginning
March 2020 which was accompanied by a spike in market volatility. We analyse
two sampling intervals: daily, end-of-day price changes (for which our data
availability covers the whole of 2019 and 2020) and 30-minute intervals (for
which our data is limited to the period between January 2019 up to and
including July 2020).

TAQ Database: We use intraday data from the TAQ (Trade and Quote) database. To
account for microstructure related issues, such as the bid-ask-bounce or
infrequent trading, we rely upon quoted prices that we re-sample to a
30-minute frequency.111When we talk about quotes, we refer to the National
Best Bid and Offer (NBBO) where the national best bid (offer) is the best
available quoted bid (offer) price across all U.S. exchanges. See [28] for an
overview. For that, we first remove all crossed quotes, i.e., all quotes where
the bid price exceeds the ask, require the bid-ask-spread to be below 5$, and
finally use the last valid available quote within every 30-minute
interval.222We forward fill quotes if there is no valid entry for a given time
interval. However, this does almost never happen for very liquid stocks such
as those chosen in this paper. We further account for dividend payments and
stock splits, which mechanically influence stock prices, and create a
performance price index using quoted mid-prices.

CRSP: We also consider lower-frequency, daily, data from the Center of
Research in Security Prices (CRSP), which is one of the most used database in
economics and finance. We again calculate a performance price index for each
stock using the daily holding period return provided by CRSP. Note that, while
we use quoted mid-prices for the 30-minute high-frequency data, CRSP uses
trade prices to calculate the holding period return. However, as the trading
volume has increased considerably over the last decade, this should not be an
issue [29].

Company | Business Sector | Ticker
---|---|---
Apple | Technology | AAPL
Citigroup | Banking | C
Walt Disney Co. | Media | DIS
Evergy Inc. | Energy | EVRG
General Electrics | Industry | GE
Pfizer | Pharmaceutics | PFE
Walmart Inc. | Retail | WMT
Table 1: The companies whose data has been analysed in our article.

##  III Theoretical Background and Model

The standard stochastic differential equation to describe a stock price $P$ is
the Geometric Brownian Motion given by

$\odv{P}{t}=\mu P+\sigma P\epsilon$ (1)

with standard Gaussian noise
$\epsilon=\epsilon(t)\overset{iid}{\sim}\mathcal{N}(0,1)$, constant drift
$\mu$ (typically $\mu>0$) and volatility $\sigma$. As pointed out in [20], the
physical interpretation as a particle’s trajectory $P$ in a potential $V(P)$
transforms equation (1) to

$\odv{P}{t}=-\odv{V}{P}(P)+\sigma P\epsilon_{t}\quad\textrm{with}\quad
V(P)=-\frac{\mu}{2}P^{2}$ (2)

and an arbitrary constant $C$ (set to zero for simplicity). However, analysing
this potential $V$ in terms of its linear stability (cf. e.g. [30]) leads to
the problematic result that the only fixed point in the data with
$\odv{V}{P}(P_{0})=0$, namely $P_{0}=0$, is an unstable fixed point for
$\mu>0$. Without a stable fixed point, trajectories are expected to diverge
away from $P_{0}=0$ towards infinity. As this is - at least for limited time
scales - a highly unrealistic model, the authors of [20, 18] have suggested
higher order polynomials in the potential $V$ of (2). From the assumption that
the rate of capital injection by investors should depend on the current market
capitalisation, they derive a quartic potential

$\displaystyle V(P)$ $\displaystyle=$
$\displaystyle-P\left(\frac{\alpha_{1}}{2}P+\frac{\alpha_{2}}{3}P^{2}+\frac{\alpha_{3}}{4}P^{3}\right)\textmd{
with }$ $\displaystyle-\odv{V}{P}(P)$ $\displaystyle=$
$\displaystyle\alpha_{1}P+\alpha_{2}P^{2}+\alpha_{3}P^{3},$ (3)

i.e. the drift term $-\odv{V}{P}(P)$ is a polynomial of order $q=3$. With a
suitable choice of parameters $\alpha$, this potential can adopt the shape of
a double-well potential with stable fixed points at $P_{0}=0$ and $P_{1}>0$,
thereby predicting both the presence of a bankruptcy state at the stable fixed
point $P_{0}=0$ and an additional stable state with nonzero price $P_{1}>0$.
In [20], however, major constraints to the parameters during the estimation
process were necessary to achieve this.

###  III.1 Numerical Implementation

If a time series $(s_{n}(t_{n}))_{n}$ with observations $s_{n}$ at time
$t_{n}$ has been recorded, a maximum likelihood approach can be used to
estimate the most likely parameters $\sigma$ and $\alpha_{i}$ such that a
stochastic differential equation according to equations (2) and (III) may have
produced the observed time series. For any two adjacent points $s_{n}(t_{n})$
and $s_{n+1}(t_{n+1})$ and any given parameters
$\phi=(\sigma^{2},\alpha_{i})$, the likelihood $L$ of observing the transition
from $s_{n}(t_{n})$ to $s_{n+1}$ at $t_{n+1}$ can be explicitly calculated as

$\displaystyle L(s_{n+1}|t_{n+1},s_{n},t_{n},\phi)=\left(2\pi\left(\sigma
s_{n}\sqrt{t_{n+1}-t_{n}}\right)^{2}\right)^{-\frac{1}{2}}$ (4)
$\displaystyle\cdot\exp{\left(-\frac{\left(s_{n+1}-\left(s_{n}+\left(-\odv{V}{s}(s)\right)(t_{n+1}-t_{n})\right)\right)^{2}}{2\left(\sigma
s_{n}\sqrt{t_{n+1}-t_{n}}\right)^{2}}\right)}$

or as the log likelihood

$\displaystyle\mathcal{L}(s_{n+1}|t_{n+1},s_{n},t_{n},\phi)=\log
L(s_{n+1}|t_{n+1},s_{n},t_{n},\phi)$ (5)
$\displaystyle=-\frac{1}{2}\log\left(2\pi\left(\sigma
s_{n}\sqrt{t_{n+1}-t_{n}}\right)^{2}\right)$
$\displaystyle-\frac{\left(s_{n+1}-\left(s_{n}+\left(-\odv{V}{s}(s)\right)(t_{n+1}-t_{n})\right)\right)^{2}}{2\left(\sigma
s_{n}\sqrt{t_{n+1}-t_{n}}\right)^{2}}.$

Because we assume Markovian dynamics, the complete log likelihood for the full
observed time series is then simply the sum over the stepwise log likelihoods

$\mathcal{L}\left((s_{n})_{n}|(t_{n})_{n},\phi\right)=\sum_{i=0}^{n-1}\mathcal{L}(s_{i+1}|t_{i+1},s_{i},t_{i},\phi).$
(6)

For given observations $(s_{i},t_{i})$, the likelihood $\mathcal{L}$ can be
maximised by varying the parameters $\phi$ to estimate the optimal parameters
$\phi^{*}$.

According to Bayes’ Theorem and Bayesian Statsitics [31], the likelihood of
observing the measured data conditional on some parameter values
$L\left(\left(s_{n}(t_{n})\right)_{n}\middle|\phi\right)$ is combined with an
a-priori distribution $f_{prior}(\phi)$ to calculate a posterior distribution
of the parameters given the observed data:

$f_{post}(\phi|(s_{n}(t_{n}))_{n})\sim
f_{prior}(\phi)L\left(\left(s_{n}(t_{n})\right)_{n}\middle|\phi\right).$ (7)

For an uninformed flat prior, this transformation is mathematically trivial,
but allows us to calculate $f_{post}$ as a probability density of the
parameters $\phi$ conditional on the observed data. Hence, the distribution of
the parameters $\phi$ can be explored via Markov chain Monte Carlo (MCMC)
methods (e.g. [32]) by drawing samples $(\phi^{(j)})_{j}$ from the posterior
distribution as implemented in the Python package emcee [33]. MCMC can uncover
correlations between different parameters and also explore local maxima of the
probability density. It therefore gives a more complete view of the underlying
distribution than summary statistics like e.g. the mean or standard deviation.
In particular, we will use the sampled parameters $(\phi^{(j)})_{j}$ to
construct an ensemble of potentials $V(P)$ and evaluate whether their shapes
are roughly consistent with each other.

###  III.2 Synthetic Data

To test our method, synthetic time series $(s_{n})_{n}$ are simulated via the
Euler-Maruyama scheme [34] as

$s_{n+1}=s_{n}+\left(-\odv{V}{s}(s)\right)(t_{n+1}-t_{n})+\sigma
s_{n}\sqrt{t_{n+1}-t_{n}}\epsilon_{n}$ (8)

with $\epsilon_{n}\overset{iid}{\sim}\mathcal{N}(0,1)$ for any parameters
$\alpha_{i}$ for the potential in (III). The following paragraphs discuss how
well our model can then identify the underlying dynamcis and parameters from
the observed data $(s_{n},t_{n})$. Note that for the synthetic data, non-
equidistant time steps have been used.

####  III.2.1 Estimating the Correct Order

Generalising the potential from (III) to a potential with arbitrary polynomial
order $q$ leads to

$V(P)=-P\sum_{i=1}^{q}\frac{\alpha_{i}}{i+1}P^{i}.$ (9)

For given $q$, random parameter values $\alpha_{i}$ and a random noise level
$\sigma$ are sampled and the resulting time series is simulated with these
parameters. If the sampled parameters result in numerical errors (i.e. an
infinitely diverging time series), the time series is discarded from the
ensemble. This is done repeatedly until the ensemble includes 100 time series
with a length of 1000 time steps each. For those synthetic time series, the
best order is estimated via the Akaike information criterion $AIC$ [35] given
by

$AIC=-2\mathcal{L}_{\max}+2(q+1)$ (10)

where $q+1$ is the total number of model parameters ($q$ monomials’ prefactors
$\alpha_{i}$ and $\sigma$). Hence, $q$ is varied, the maximum likelihood
$\mathcal{L}_{\max}$ for the chosen $q$ is estimated and the resulting $AIC$
is calculated. The model with the lowest $AIC$ is chosen as the best model.
The results shown in figure 1 show that for polynomial orders $q=1$ (Geometric
Brownian Motion), $q=2$ and $q=3$ (Halperin’s suggestion as in equation
(III)), the correct order is usually identified as such. An even higher order
$q=4$ shows very unreliable results, but will be included for completeness for
the further data analysis.

Figure 1: For each true polynomial order $q$, 100 trajectories are randomly
sampled and then their best order is estimated. The histograms show that the
method successfully estimates trajectories with order $q=1,2,3$, but struggles
with the higher order $q=4$.

####  III.2.2 Estimating the Parameters

Instead of sampling repeated trajectories with different parameters, now for
order $q=3$, the parameters
$\phi=(\sigma^{2},\alpha_{1},\alpha_{2},\alpha_{3})$ are kept constant as
$\phi=(0.05,2,-1,0.01)$. 100 time series with 1000 time steps are sampled and
their parameters are estimated by fitting a model with $q=3$. Despite the
constant parameters, the randomness of the $\epsilon_{n}$ in equation (8)
nevertheless ensures that the time series are different from each other. The
histograms of the estimated parameters, their means and standard deviations
are shown together with the true parameter values in figure 2. Note that the
true parameter value is always within the one-standard-deviation-interval
around the mean and that the parameter $\alpha_{3}$ has a distribution
virtually indistinguishable from that of a parameter with mean zero: The
parameter estimation correctly shows that $\alpha_{3}$ is so low that it is a
superfluous parameter for model inference. Note that if the same data is
estimated by a model with $q=2$, the results are fairly consistent with the
depicted histograms. However, fitting the data to a model with $q=4$ results
in the one-standard-deviation-interval of $\alpha_{2}$ also containing the
value 0, which is a consequence of overfitting the model.

Figure 2: Parameter estimations for 100 trajectories with the same true
parameters $\phi=(\sigma^{2},\alpha_{i})$ as given by the solid lines. The
interval of Mean $\pm$ Standard Deviation of the estimated ensembles always
includes the true parameters.

##  IV Results

While [20] analyses time periods of one year to estimate parameters, we
believe that because of the assumption of constant volatility in equation (4),
it is prudent to restrict the date to shorter intervals of one trading month.
Hence, we divide the given data into non-overlapping monthly intervals and
estimate the polynomial order $q$ of the underlying stochastic differential
equation via the $AIC$. Note that the time difference between each observation
is taken as a constant interval of 1 time step in trading days or 30-minute-
steps, respectively, including the overnight return.

###  IV.1 Polynomial Orders

The distributions of the estimated polynomial orders $q$ are shown in figure
3: On both time scales and in all market periods, the order $q=2$ is the most
frequently estimated order with the GBM model at $q=1$ being the second most
frequent estimation. The suggestion $q=3$ from [20] as well as the even
higher-order $q=4$ are only rarely estimated as the most accurate model.
Interestingly, calm and turbulent periods (as defined in section II) show
essentially identical distributions, whereas the order $q=4$ seems to be a bit
more frequently estimated for the shorter time scale of 30 minutes than for
the daily data. Table 2 shows that there is high consistency between the
estimated orders for both time intervals for orders $q=1$ and $q=2$, but
increasing disagreements for orders $q=3$ and $q=4$.
Overall, this suggests that a polynomial order of $q=1$ and $q=2$ can be a
reasonable modelling assumption for the time series data and that the
identification of these two orders is consistent for the two sampling
intervals under consideration, whereas the general market state does not seem
to influence our results.

Figure 3: Estimated polynomial orders $q$ of the monthly real-world price time series of the companies in table 1 for the daily (left) and 30-minute-intervals (right). Differences between the calm (up to February 2020) and turbulent period (starting in March 2020) are only small.  | Optimal Order
---
(30 Min)
Opt. Order (Daily)
| 1 | 2 | 3 | 4
1 | 24 | 6 | 3 | 1
2 | 6 | 41 | 10 | 3
3 | 1 | 10 | 8 | 1
4 | 7 | 7 | 3 | 2
Table 2: Comparison of the estimated orders $q$ for the same company and month
with the price time series in daily and 30-minute-intervals.

###  IV.2 Potentials

For the optimal polynomial order $q$, we then sampled the parameters
$\phi=(\sigma^{2},\alpha_{1},\dots,\alpha_{q})$ from their posterior
probability distribution to get an ensemble of parameters $(\phi^{(j)})_{j}$.
From that, we calculate the corresponding ensemble of potentials $V(P)$
according to (9) and plot them, their pointwise centred 68% and 95% credible
intervals (CIs) and the potential corresponding to the maximum likelihood
estimation. The zero horizontal is shown in these plots as the y-axis position
of the potential at $P=0$ to indicate where the potential is above or below
the potential energy at the zero price (and if the price ”particle” would
therefore prefer or not prefer to be at the potential level of $P=0$). A
couple of generic features can be observed for these potentials and do not
depend on the chosen sampling rate:

####  IV.2.1 Order 2

As the order $q=2$ is the most frequently identified polynomial order
according to the results in 3, it is quite insightful to focus on the
associated potentials. They virtually always look like the potential depicted
in figure 4 and show a potential well as a pronounced minimum. Close to this
minimum, the 68% CI is usually also below 0 and sometimes (as depicted in
figure 4) even the 95% CI. The MCMC-sampled potential ensembles thus support
the existence of a potential well as they clearly show the potential well for
a large majority of trajectories. Following the interpretation of the
potential wells from [20], this supports the existence of a locally stable
price within this potential minimum.

Figure 4: MCMC-sampled potentials of order $q=2$ for 30-minute-intervals for
Pfizer in September 2019. The maximum likelihood estimation (Estimation) lies
firmly within the ensemble of potentials and even the 68% credible interval
shows a clear potential well.

####  IV.2.2 Order 1: GBM

The GBM with $q=1$ is the second most frequently estimated order. As shown in
the subplots $a$ and $b$ in figure 5, the MCMC samples show two general types
of ensembles: in $a$, the maximum likelihood estimation of the potential is
always very close to 0 and the 68% CIs therefore envelope the zero horizontal.
This makes it difficult to gauge a clear direction of the potential and hence
of the movement of the price time series. Contrary to that, potential
ensembles like in $b$ have a clear direction. In $b$, the potential is
increasing for higher prices (and hence a restoring force pulls the price to
the minimum at 0), but decreasing potential ensembles can also be found for
other time intervals. That means some time intervals like $b$ show a
predominant direction of price movement, whereas others like $a$ have no
predominant direction, but rather a random movement.

Figure 5: MCMC-sampled potentials of order $q=2$ for 30-minute-intervals for
Pfizer in time intervals with $q=1$ ($a$ and $b$), $q=3$ ($c$) and $q=4$
($d$). Note that the maximum likelihood estimation lies within the $68\%$
credible interval for all orders $q$ except for $q=4$.

####  IV.2.3 Order 3

The potential with $q=3$ is the one suggested in [20] and the typical shape of
their MCMC samples are shown in subplot $c$ of figure 5. Note that some of
these potentials are also mirrored along the x-axis. Similar to the potentials
in $b$, they also show a predominant direction, but also often a bistable
saddle point. Notably, they do not show the pronounced double-well potential
predicted in [20].

####  IV.2.4 Order 4

Potentials with $q=4$ usually correspond to very wide potential wells as can
be seen in subplot $d$ of figure 5 (compare e.g. its full width at half
maximum to that of the potential in 4). In the depicted MCMC ensemble, the
maximum likelihood estimation does not lie within the 68% credible interval.
This indicates a multimodal posterior distribution and was found in
surprisingly many time intervals. Similar to the ensemble shown in subfigure
$a$, these potentials also often envelope the x-axis with their 68% CIs and
therefore show no clearly predominant direction.

##  V Conclusion

###  V.1 Summary

We use a maximum likelihood estimation to analyse price time series of stocks.
Via the $AIC$ model selection, we find that a second order polynomial for the
drift term often offers a very suitable description of the data. While the
standard GBM model with a first order polynomial is not selected as frequently
as the second order model, it still appears often enough to be considered a
valid candidate model. Higher order polynomials are rarely estimated. Sampling
the posterior density of the parameters via MCMC reveals that the second order
polynomials’ potentials show very pronounced potential wells (i.e. stable
minima) for nonzero prices which is mathematically impossible for the GBM’s
potential as pointed out in [20].

###  V.2 Discussion

Our research question is heavily inspired by [20], but differs from it in a
key factor: The model presented in [20] always has a drift polynomial of order
$q=3$ and uses the credit default swap rates (CDS) to estimate the probability
of a considered company going bankrupt. This probability is then used as a
constraint in the parameter estimation such that the jump from a potential
well with nonzero price to another potential well at price zero (i.e. the
stock collapsing) has a jump probability (Kramers’s escape rate) which is
quantified by the CDS. Thus, the work in [20] combines the price dynamics of
stochastic differential equations with the CDS data as additional constraints
to estimate a stochastic differential equation with a probability of the stock
crashing. In contrast to this, our estimation scheme uses no additional
constraints or external data, but purely the price time series. Our $q=3$
estimations (subfigure $c$ in 5) do not show the double-well potential
postulated by [20]. However, our model selection via the $AIC$ indicates that
$q=2$ is instead the most frequently observed polynomial degree and for $q=2$,
MCMC shows a clear potential well that is consistent for the whole MCMC
ensemble. In short, because we do not use additional constraints, we cannot
reproduce the double well potential with default probabilities, but instead
show via our fully unconstrained approach that the potential wells arise
naturally just from the price time series alone.
It is also noteworthy that the $AIC$ selection means that the stocks are not
described by the same polynomial degree $q$ for all time series. Instead, our
model selection implies that the potential itself is time-dependent. A
possible explanation for this is that due to external effects, an order
parameter such as the capital influx into the financial markets is changed.
Then, the underlying potential might change due to these effects and e.g.
experience a bifurcation, resulting in changing price dynamics. Whether a
bifurcation is a suitable description of the dynamics under such conditions
requires further analysis on the transitions between the different models. One
can imagine e.g. that a price initially starts as being in a stable
equilibrium with $q=2$ like in figure 4, but external news change the
potential to that of subfigure $b$ in figure 5 with $q=1$. Now, the price has
a predominant direction of movement and is no longer experiencing a restoring
force back to the equilibrium price and can therefore explore new areas of the
phase space. The market finally manages to process the news and their
implications and finally, the price reaches a new equilibrium with $q=2$.
Thus, the price at the potential minimum can be interpreted as a fair price
similar to the discussion in [17]. However, further research into the
transition between the different potentials is necessary to verify this
interpretation. Note that GBMs with potentials such as subfigure $a$ in figure
5 have essentially no predominant direction of movement and show a random walk
without a restoring force. This is a different behaviour to $q=2$ which also
does not show a predominant direction, but instead has such a restoring force
that restricts the price to the potential well.
One might have assumed that stable fixed points ($q=2$) should occur
significantly less frequently during the turbulent period because of the
overall instability of the market. But interestingly, our results do not seem
to show a difference between the calm and turbulent market period (cf. figure
3), perhaps indicating that the market can quickly adjust to such turbulent
behaviour.
Finally, it is a reassuring result that the $AIC$ selection still frequently
suggests $q=1$ (the standard GBM model) as the best polynomial order. The
standard GBM model still appears rather frequently in our data and therefore
nevertheless manages to provide a reasonably accurate model.

###  V.3 Further Research and Applications

While this article focused on the drift term like in [20], there are of course
possible extensions of the diffusion/volatility that can be taken into
account, too. Stochastic volatility and local volatility models have been
widely accepted in finance [36, 37], but other modelling possibilities exist,
too: While our article used the volatility parametrisation from the GBM in (1)
via $\sigma P\epsilon_{t}$, one can also imagine e.g. a polynomial model here
given by

$\textmd{Diffusion}(P)=\sigma\left(\sum_{j}\beta_{j}P^{j}\right)\epsilon_{t}.$
(11)

However, from the author’s experience, the maximum likelihood estimation can
become troublesome if the diffusion term has several free parameters as the
estimator can then attempt to essentially attribute the whole observed
dynamics to the diffusion. A strong regularisation might be necessary if one
wishes to expand the diffusion model.
Another model extension might be the incorporation of memory effects.
Generalised versions of the Langevin equation include non-Markovian memory
terms by e.g. an explicit memory kernel [38] or by assuming the existence of a
second hidden process that has not been observed [39]. Such a hidden component
might correspond to the traders’ knowledge or belief which certainly
influences the stock prices, but is not explicitly recorded.
Finally, it is noteworthy to point out that although we used equidistant time
intervals between the observations of the data, the model has been tested on
synthetic data with non-equidistant time intervals in section III.2. Such a
situation arises naturally in the context of tick-by-tick data which is the
highest resolution of trading data. Here, instead of sampling the price at a
high frequency, every single trade is recorded at the exact time that it
occurred. As the time between two subsequent trades can be arbitrarily short
or long, the application of a robust method without the need for equidistant
time steps might prove useful here.

###### Acknowledgements.

The authors thank Tim Kroll (WWU Münster) for valuable discussion about the
propagator and its technical implementation. Tobias Wand is financed by the
Studienstiftung des deutschen Volkes.

## References

  * Mantegna [1997] R. N. Mantegna, Degree of correlation inside a financial market, AIP Conference Proceedings 411, 197 (1997).
  * Laloux  _et al._ [1999] L. Laloux, P. Cizeau, J.-P. Bouchaud, and M. Potters, Noise dressing of financial correlation matrices, Phys. Rev. Lett. 83, 1467 (1999).
  * Münnix  _et al._ [2012] M. C. Münnix, T. Shimada, R. Schäfer, F. Leyvraz, T. H. Seligman, T. Guhr, and H. E. Stanley, Identifying states of a financial market, Scientific Reports 2 (2012).
  * Mantegna and Stanley [1995] R. N. Mantegna and H. E. Stanley, Scaling behaviour in the dynamics of an economic index, Nature 376, 46 (1995).
  * Lillo  _et al._ [2005] F. Lillo, S. Mike, and J. D. Farmer, Theory for long memory in supply and demand, Physical Review E 71, 10.1103/physreve.71.066122 (2005).
  * Sato and Kanazawa [2023] Y. Sato and K. Kanazawa, Direct quantitative evidence of the order-splitting hypothesis as the microscopic origin of long-range correlations in market order flow (2023), arXiv:2301.13505 [q-fin.TR] .
  * Mantegna and Stanley [2000] R. Mantegna and H. Stanley,  _An Introduction to Econophysics_ (Cambridge University Press, Cambridge, 2000).
  * Da Cunha [2021] C. Da Cunha,  _Introduction to Econophysics: Contemporary Approaches with Python Simulations_ (CRC Press, 2021).
  * Bachelier [1900] L. Bachelier, Théorie de la spéculation, Annales scientifiques de l’École Normale Supérieure 3e série, 17, 21 (1900).
  * Black and Scholes [1973] F. Black and M. Scholes, The pricing of options and corporate liabilities, Journal of Political Economy 81, 637 (1973).
  * Merton [1973] R. C. Merton, Theory of rational option pricing, The Bell Journal of Economics and Management Science 4, 141 (1973).
  * Renner  _et al._ [2001] C. Renner, J. Peinke, and R. Friedrich, Evidence of markov properties of high frequency exchange rate data, Physica A: Statistical Mechanics and its Applications 298, 499 (2001).
  * Tice and Webber [1997] J. Tice and N. Webber, A nonlinear model of the term structure of interest rates, Mathematical Finance 7, 177 (1997), https://onlinelibrary.wiley.com/doi/pdf/10.1111/1467-9965.00030 .
  * Bouchaud and Cont [1998] J.-P. Bouchaud and R. Cont, A langevin approach to stock market fluctuations and crashes, Physics of Condensed Matter 6, 543 (1998).
  * Kanazawa  _et al._ [2018] K. Kanazawa, T. Sueshige, H. Takayasu, and M. Takayasu, Kinetic theory for financial brownian motion from microscopic dynamics, Phys. Rev. E 98, 052317 (2018).
  * Lee  _et al._ [2002] M. Lee, A. Oba, and H. Takayasu, Parameter estimation of a generalized langevin equation of market price (Springer, 2002) pp. 260–270.
  * Garcia  _et al._ [2020] M. Garcia, A. Machado Pereira, J. Acebal, and A. Bosco de Magalhães, Forecast model for financial time series: An approach based on harmonic oscillators, Physica A: Statistical Mechanics and its Applications 549, 124365 (2020).
  * Halperin [2020] I. Halperin, The inverted parabola world of classical quantitative finance: Non-equilibrium and non-perturbative finance perspective (2020), arXiv:2008.03623 [q-fin.GN] .
  * Dixon  _et al._ [2020] M. F. Dixon, I. Halperin, and P. Bilokon,  _Machine Learning in Finance_ (Springer International Publishing, 2020).
  * Halperin and Dixon [2019] I. Halperin and M. Dixon, “quantum equilibrium-disequilibrium”: Asset price dynamics, symmetry breaking, and defaults as dissipative instantons, Physica A: Statistical Mechanics and its Applications 537, 122187 (2019).
  * Siegert  _et al._ [1998] S. Siegert, R. Friedrich, and J. Peinke, Analysis of data sets of stochastic systems, Physics Letters A 243, 275 (1998).
  * Friedrich  _et al._ [2000a] R. Friedrich, S. Siegert, J. Peinke, S. Lück, M. Siefert, M. Lindemann, J. Raethjen, G. Deuschl, and G. Pfister, Extracting model equations from experimental data, Physics Letters A 271, 217 (2000a).
  * Friedrich  _et al._ [2000b] R. Friedrich, J. Peinke, and C. Renner, How to quantify deterministic and random influences on the statistics of the foreign exchange market, Phys. Rev. Lett. 84, 5224 (2000b).
  * Kleinhans and Friedrich [2007] D. Kleinhans and R. Friedrich, Maximum likelihood estimation of drift and diffusion functions, Physics Letters A 368, 194 (2007).
  * Brunton  _et al._ [2016] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 113, 3932 (2016).
  * Fasel  _et al._ [2022] U. Fasel, J. N. Kutz, B. W. Brunton, and S. L. Brunton, Ensemble-SINDy: Robust sparse model discovery in the low-data, high-noise limit, with active learning and control, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 478 (2022).
  * Heßler and Kamps [2021] M. Heßler and O. Kamps, Bayesian on-line anticipation of critical transitions, New Journal of Physics (2021).
  * Holden and Jacobsen [2014] C. W. Holden and S. Jacobsen, Liquidity measurement problems in fast, competitive markets: Expensive and cheap solutions, The Journal of Finance 69, 1747 (2014).
  * Wiedemann [2022] T. Wiedemann, Quantifying the beta estimation bias and its implications for empirical asset pricing, Available at SSRN 4179365 (2022).
  * Strogatz [2000] S. Strogatz,  _Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering_ , Studies in nonlinearity (Westview, 2000).
  * Von der Linden  _et al._ [2014] W. Von der Linden, V. Dose, and U. Von Toussaint,  _Bayesian probability theory: applications in the physical sciences_ (Cambridge University Press, 2014).
  * Hastings [1970] W. K. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika 57, 97 (1970).
  * Foreman-Mackey  _et al._ [2013] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: The MCMC hammer, Publications of the Astronomical Society of the Pacific 125, 306 (2013).
  * Kloeden and Platen [2013] P. E. Kloeden and E. Platen,  _Numerical solution of stochastic differential equations_ , 1992nd ed., Stochastic Modelling and Applied Probability (Springer, Berlin, Germany, 2013).
  * Akaike [1998] H. Akaike, Information theory and an extension of the maximum likelihood principle, in  _Selected Papers of Hirotugu Akaike_, edited by E. Parzen, K. Tanabe, and G. Kitagawa (Springer New York, New York, NY, 1998) pp. 199–213.
  * Heston [1993] S. L. Heston, A closed-form solution for options with stochastic volatility with applications to bond and currency options, The review of financial studies 6, 327 (1993).
  * Dupire  _et al._ [1994] B. Dupire _et al._ , Pricing with a smile, Risk 7, 18 (1994).
  * Willers and Kamps [2021a] C. Willers and O. Kamps, Efficient bayesian estimation of the generalized Langevin equation from data (2021a), arXiv:2107.04560 [physics.data-an] .
  * Willers and Kamps [2021b] C. Willers and O. Kamps, Non-parametric estimation of a langevin model driven by correlated noise, The European Physical Journal B 94 (2021b).


# Poisson-Gamma Dynamical Systems with Non-Stationary Transition Dynamics

Jiahao Wang1 Sikun Yang2 Heinz Koeppl3 Xiuzhen Cheng1 Pengfei Hu1 Guoming
Zhang1 1School of Computer Science and Technology, Shandong University
2School of Computing and Information Technology, Great Bay University
3Department of Electrical Engineering and Information Technology, Technische
Universität Darmstadt

###### Abstract

Bayesian methodologies for handling count-valued time series have gained
prominence due to their ability to infer interpretable latent structures and
to estimate uncertainties, and thus are especially suitable for dealing with
_noisy_ and _incomplete_ count data. Among these Bayesian models, Poisson-
Gamma Dynamical Systems (PGDSs) are proven to be effective in capturing the
evolving dynamics underlying observed count sequences. However, the state-of-
the-art PGDS still falls short in capturing the _time-varying_ transition
dynamics that are commonly observed in real-world count time series. To
mitigate this limitation, a non-stationary PGDS is proposed to allow the
underlying transition matrices to evolve over time, and the evolving
transition matrices are modeled by sophisticatedly-designed Dirichlet Markov
chains. Leveraging Dirichlet-Multinomial-Beta data augmentation techniques, a
fully-conjugate and efficient Gibbs sampler is developed to perform posterior
simulation. Experiments show that, in comparison with related models, the
proposed non-stationary PGDS achieves improved predictive performance due to
its capacity to learn non-stationary dependency structure captured by the
time-evolving transition matrices.

##  1 Introduction

In recent years, there has been an increasing interest in modeling count time
series. For instance, some previous works Blei and Lafferty (2006); Wang and
McCallum (2006); Jähnichen et al. (2018) are concerned with how to learn the
evolving topics behind text corpus (frequencies of words) over time. Some
works Sheldon and Dietterich (2011); Raymer et al. (2013); Wilson (2017);
Wanner (2021) try to predict global immigrant trends underlying international
population movements. Count time series are often _overdispersed_ , _sparse_ ,
_high-dimensional_ , and thus can not be well modeled by widely used dynamic
models such as linear dynamical systems Kalman (1960); Ghahramani and Roweis
(1998). Recently, many works Zhou and Carin (2012, 2015); Schein et al. (2015,
2016b, 2016a); Acharya et al. (2015); Schein et al. (2019) prefer to choose
distributions of the gamma-Poisson family to build their hierarchical Bayesian
models. In particular, these models enjoy strong explainability and can
estimate uncertainty especially when the observations are _noisy_ and
_incomplete_. Among these works, Poisson-Gamma Dynamical Systems (PGDSs)
Schein et al. (2016b) received a lot of attention because PGDS can learn how
the latent dimensions excite each other to capture complicated dynamics in
observed count series. For instance, a very inspiring research paper may
motivate other researchers to publish papers on related topics Chang and Blei
(2009). The outbreak of COVID-19 in one state, may lead to the rapid rising of
COVID-19 cases in the nearby states and vice versa Unwin et al. (2020). In
particular, PGDS can be efficiently learned with a tractable Gibbs sampling
scheme via Poisson-Logarithmic data augmentation and marginalization technique
Zhou and Carin (2015). Due to its strong flexibility, PGDS achieves better
performance in predicting missing entities and future observations, compared
with related models Ghahramani and Roweis (1998); Acharya et al. (2015).

Despite these advantages, PGDS still can not capture the time-varying
transition dynamics underlying observed count sequences, which are commonly
observed in real-world scenarios Winkelmann (2008). For instance, during the
initial stage of the COVID-19 pandemic, the worldwide counts of infectious
patients were significantly affected by various local policies, government
interventions, and emergent events Grossman et al. (2020); Team (2021); Feng
et al. (2021). The cross transition dynamics among the different monitoring
areas were also evolving as the corresponding policies and interventions
changed over time. Hence, PGDS unavoidably makes a certain amount of
approximation error in capturing the aforementioned non-stationary count time
series, using a _time-invariant_ transition kernel.

To mitigate this limitation, Non-Stationary Poisson-Gamma Dynamical Systems
(NS-PGDSs), a novel kind of Poisson-gamma dynamical systems with non-
stationary transition dynamics are developed. More specifically, NS-PGDS
captures the evolving transition dynamics by sophisticatedly designed
Dirichlet Markov chains. Via the Dirichlet-Multinomial-Beta data augmentation
strategy, the Non-Stationary Poisson-Gamma Dynamical Systems can be inferred
with a conjugate-yet-efficient Gibbs sampler. Our contributions are summarized
as follows:

  * •

We propose a Non-Stationary Poisson-Gamma Dynamical System (NS-PGDS), a novel
Poisson-gamma dynamical system with time-evolving transition matrices that can
well capture non-stationary transition dynamics underlying observed count
series.

  * •

Three sophisticated Dirichlet Markov chains are dedicated to improving the
flexibility and expressiveness of NS-PGDSs, for capturing the complex
transition dynamics behind sequential count data.

  * •

Fully-conjugate-yet-efficient Gibbs samplers are developed via Dirichlet-
Multinomial-Beta augmentation techniques to perform posterior simulation for
the proposed Dirichlet Markov chains.

  * •

Extensive experiments are conducted on four real-world datasets, to evaluate
the performance of the proposed NS-PGDS in predicting missing and future
unseen observations. We also provide exploratory analysis to demonstrate the
explainable latent structure inferred by the proposed NS-PGDS.

##  2 Preliminaries

Let
$\bm{y}^{\left(t\right)}=\left[y_{1}^{(t)},\cdots,y_{V}^{(t)}\right]^{\mathrm{T}}\in\mathbb{N}_{0}^{V}$
be a vector of nonnegative count valued observations at time $t$. To capture
the latent dynamics underlying count sequences, some previous works Han et al.
(2014); Kalantari and Zhou (2020) model the observations as

$\bm{y}^{(t)}=p\left(\bm{z}^{(t)}\right),\ \
\bm{z}^{(t)}=f^{-1}\left(\bm{x}^{(t)}\right),$

where $p\left(\cdot\right)$ is the observation likelihood function, and
$f\left(\cdot\right)$ is a link function that maps the parameters of
observation component to continuous-valued latent variables
$\bm{x}^{(t)}\in\mathbb{R}^{K}$. The latent factor $\bm{x}^{(t)}$ evolves over
time according to a linear dynamical system given by
$\bm{x}^{(t)}\sim\mathcal{N}\left(\bm{A}\bm{x}^{(t-1)},\bm{\Lambda}^{-1}\right)$,
where $\bm{A}$ is the state transition matrix of size $K\times K$, and
$\bm{\Lambda}=\mathrm{diag}\left(\lambda_{1},\cdots,\lambda_{K}\right)$ is the
inverse covariance matrix with $\lambda_{k}^{-1}$ determining the variance of
$k$-th latent dimension. Han et al. Han et al. (2014) adopted the Extended
Rank likelihood function to model count observations using LDS with time
complexity $\mathcal{O}\left(\left(K+V\right)^{3}\right)$, which prevents it
from practical applications for analyzing large-scale count data.

Recently, Acharya et al. Acharya et al. (2015) and Schein et al. Schein et al.
(2016b, 2019) developed Poisson-gamma family models for sequential count
observations. Gamma Process Dynamic Poisson Factor Analysis (GP-DPFA) Acharya
et al. (2015) models count data as
$y_{v}^{(t)}\sim\mathrm{Pois}\left(\sum_{k=1}^{K}\lambda_{k}\phi_{vk}\theta_{k}^{(t)}\right)$,
where $\theta_{k}^{(t)}$ represents the strength of $k$-th latent factor at
time $t$, and $\phi_{vk}$ captures the involvement degree of $k$-th factor to
$v$-th observed dimension. To ensure the model identifiability, we can impose
a restriction as $\sum_{v}\phi_{vk}=1$, and thus place a Dirichlet prior over
$\bm{\phi_{k}}=\left[\phi_{1k},\cdots,\phi_{Vk}\right]^{T}$ as
$\bm{\phi_{k}}\sim\mathrm{Dir}\left(\epsilon_{0},\cdots,\epsilon_{0}\right)$.
To capture the underlying dynamics, the latent factor $\theta_{k}^{(t)}$
evolves over time according to a gamma Markov chain as
$\theta_{k}^{(t)}\sim\mathrm{Gam}\left(\theta_{k}^{(t-1)},c_{t}\right)$, where
$c_{t}$ is the rate parameter of the gamma distribution to control the
variance of gamma Markov chains. Although GP-DPFA can well fit one-dimensional
count sequences, it fails to learn how the latent dimensions interact with
each other.

To address this concern, Schein et al. Schein et al. (2016b) developed
Poisson-gamma dynamical systems to capture the underlying transition dynamics.
In particular, $\theta_{k}^{(t)}$ evolves over time as
$\theta_{k}^{(t)}\sim\mathrm{Gam}\left(\tau_{0}\sum_{k_{2}=1}^{K}\pi_{kk_{2}}\theta_{k_{2}}^{(t-1)},\tau_{0}\right)$,
where $\pi_{kk_{2}}$ represents how $k_{2}$-th latent factor excites the
$k$-th latent factor at next time step, and $\sum_{k=1}^{K}\pi_{kk_{2}}=1$.

##  3 Non-Stationary Poisson-Gamma Dynamical Systems

Real-world count time sequences are often _non-stationary_ because the
external interventional environments are always changing over time. The
stationary PGDS with a time-invariant transition kernel fails to capture such
time-varying transition dynamics. For instance, the transition dynamics behind
COVID-19 infectious processes are time-varying, and highly affected by various
interventional policies. Hence, to mitigate this limitation, we model the
count sequences as

$y_{v}^{(t)}\sim\mathrm{Pois}\left(\delta^{(t)}\textstyle\sum_{k=1}^{K}\phi_{vk}\theta_{k}^{(t)}\right),$

in which, the latent factors are specified by

$\theta_{k}^{(t)}\sim\mathrm{Gam}\left(\tau_{0}\textstyle\sum_{k_{2}=1}^{K}\pi_{kk_{2}}^{(t-1)}\theta_{k_{2}}^{(t-1)},\tau_{0}\right),$
(1)

where the multiplicative term
$\delta^{(t)}\sim\mathrm{Gam}\left(\epsilon_{0},\epsilon_{0}\right)$ and the
transition matrices are time-varying as
$\bm{\Pi}^{(t)}\equiv\left[\pi_{kk_{2}}^{(t)}\right]_{k,k_{2}=1}^{K}$.

Figure 1: An example illustrates the Poisson-gamma dynamical systems with non-
stationary transition kernels. The three gamma dynamic processes independently
evolve over time during the ($i-1$)-th interval. During $i$-th interval,
$\theta_{1}^{(t)}$ and $\theta_{2}^{(t)}$ gradually starts to interact with
each other while $\theta_{3}^{(t)}$ remains independent to the other two
dimensions. During ($i+1$)-th interval all the three latent components start
to interact with each other. Figure 2: Graphical representation of the NS-
PGDS. The time interval is divided into equally-spaced sub-intervals. Each
sub-interval contains $M$ time steps. The transition dynamics is stationary
within a sub-interval. In particular, the transition matrices evolve over sub-
intervals via Dirichlet Markov processes while latent factors evolve over time
steps via Eq.(1).

As shown in Figure 1, to model the time-varying transition dynamics, we assume
the whole time interval can be divided into $I$ equally-spaced sub-intervals.
The transition kernel behind complicated dynamic counts is assumed to be
_static_ within each sub-interval, while evolving over sub-intervals, to
capture non-stationary behaviours. In another word, the proposed model allows
the latent factors to evolve over time steps while the transition matrices
change over sub-intervals but assumed to be stationary within each sub-
interval, as shown in Figure 2. In particular, we let each sub-interval
contains $M$ time steps, and the $i$-th interval contains time steps
$\left\\{t\mid t=\left(i-1\right)M+1,\cdots,iM\right\\}$. We define
$i\left(t\right)$ as the function that maps time step $t$ to its corresponding
sub-interval.

Figure 3: Diagrams of the proposed Dirichlet Markov constructions. (a) is the
Dir-Dir construction. (b) is the Dir-Gam-Dir construction which takes mutation
into account. (c) illustrates the PR-Gam-Dir construction which adopts Poisson
randomized gamma distribution and can be equivalently represented as Eq.(6a)
and Eq.(6b).

Dirichlet-Dirichlet Markov processes. To capture how the underlying transition
kernel smoothly evolves over sub-intervals, we first propose the Dirichlet-
Dirichlet (Dir-Dir) Markov chain structure as

$\bm{\pi}_{k}^{(i)}\mid\bm{\pi}_{k}^{(i-1)}\sim\mathrm{Dir}\left(\eta
K\pi_{1k}^{(i-1)},\cdots,\eta K\pi_{Kk}^{(i-1)}\right),$ (2)

where $\bm{\pi}_{k}^{(i)}$ represents the $k$-th column of $\bm{\Pi}^{(i)}$,
and the prior of the scaling parameter $\eta$ is given by
$\eta\sim\mathrm{Gam}\left(e_{0},f_{0}\right)$.

The initial states are defined as
$\theta_{k}^{(1)}\sim\mathrm{Gam}\left(\tau_{0}\nu_{k},\tau_{0}\right)$. The
prior for the transition kernel of the first sub-interval is given by
$\bm{\pi}_{k}^{(1)}\sim\mathrm{Dir}\left(\nu_{1}\nu_{k},\cdots,\xi\nu_{k},\cdots,\nu_{K}\nu_{k}\right)$,
where
$\nu_{k}\sim\mathrm{Gam}\left(\frac{\gamma_{0}}{K},\beta\right)\;\mathrm{and}\;\xi,\beta\sim\mathrm{Gam}\left(\epsilon_{0},\epsilon_{0}\right)$.
Note that the expectation and variance of the transition kernel at $i$-th sub-
interval can be calculated as

$\displaystyle\mathsf{E}\left[\bm{\pi}_{k}^{(i)}\mid\bm{\pi}_{k}^{(i-1)}\right]$
$\displaystyle=\bm{\pi}_{k}^{(i-1)},$
$\displaystyle\mathsf{Var}\left[\pi_{k_{1}k}^{(i)}\mid\bm{\pi}_{k}^{(i-1)}\right]$
$\displaystyle=\frac{\pi_{k_{1}k}^{(i-1)}\left(1-\pi_{k_{1}k}^{(i-1)}\right)}{\eta
K+1},$

respectively. The transition dynamics of $i$-th sub-interval inherits the
information of the previous sub-interval, and also adapts to the data observed
in the current sub-interval. The scaling parameter $\eta$ controls the
variance of the transition matrices.

The prior specification defined in Eq.(2) by rescaling the transition matrix
at the previous sub-interval allows the transition dynamics to change
smoothly, and thus might be insufficient to capture the rapid changes observed
in complicated dynamics. To further improve the flexibility of the transition
structure, two modified Dirichlet Markov chains are studied to capture the
correlation structure between the dimensions of the transition matrices over
time.

Dirichlet-Gamma-Dirichlet Markov processes. We first introduce the Dirichlet-
Gamma-Dirichlet (Dir-Gam-Dir) Markov chain to model the evolving transition
matrices as

$\displaystyle\bm{\pi}_{k}^{(i)}$
$\displaystyle\sim\mathrm{Dir}\left(\alpha_{1k}^{(i)},\cdots,\alpha_{Kk}^{(i)}\right),$
$\displaystyle\alpha_{k_{1}k}^{(i)}$
$\displaystyle\sim\mathrm{Gam}\left(\gamma_{k}^{(i-1)}\textstyle\sum_{k_{2}=1}^{K}\psi^{(i-1)}_{kk_{1}k_{2}}\pi_{k_{2}k}^{(i-1)},c_{k}^{(i)}\right),$
(3)

where we use $\psi_{kk_{1}k_{2}}^{(i-1)}$ to capture the mutation between two
consecutive sub-intervals, and its prior is given by

$\displaystyle\left(\psi_{k1k_{2}}^{(i-1)},\cdots,\psi_{kKk_{2}}^{(i-1)}\right)$
$\displaystyle\sim\mathrm{Dir}\left(\epsilon_{0},\cdots,\epsilon_{0}\right),$
$\displaystyle\gamma_{k}^{(i)},c_{k}^{(i)}$
$\displaystyle\sim\mathrm{Gam}\left(\epsilon_{0},\epsilon_{0}\right).$

Compared with the construction defined by Eq.(2), the expectation of
Dirichlet-Gamma-Dirichlet Markov chain is

$\mathsf{E}\left[\bm{\pi}_{k}^{(i)}\mid\bm{\pi}_{k}^{(i-1)}\right]=\bm{\Psi}_{k}^{(i-1)}\bm{\pi}_{k}^{(i-1)}.$

This construction takes interactions among components of columns into account.
Hence it will dramatically improve the flexibility of our model and better fit
more complicated dynamics, compared with Dir-Dir Markov chains that only yield
smoothing transition dynamics.

Poisson-randomized-gamma-Dirichlet Markov processes. To further improve the
expressiveness of the time-varying transition kernel, and to induce sparse
patterns, we further define the Poisson-randomized-gamma-Dirichlet (PR-Gam-
Dir) Markov chain as

$\displaystyle\bm{\pi}_{k}^{(i)}$
$\displaystyle\sim\mathrm{Dir}\left(\alpha_{1k}^{(i)},\cdots,\alpha_{Kk}^{(i)}\right),$
$\displaystyle\alpha_{k_{1}k}^{(i)}$
$\displaystyle\sim\mathrm{RG1}\left(\epsilon^{\alpha},\gamma_{k}^{(i-1)}\textstyle\sum_{k2=1}^{K}\psi^{(i-1)}_{kk_{1}k_{2}}\pi_{k_{2}k}^{(i-1)},c_{k}^{(i)}\right),$
(4)

where $\mathrm{RG1}\left(\cdot\right)$ denotes the randomized gamma
distribution of the first type Yuan and Kalbfleisch (2000). Similarly, for
$\psi_{kk_{1}k_{2}}^{(i-1)}$, $\gamma_{k}^{(i)}$, and $c_{k}^{(i)}$, the
priors are given by

$\displaystyle\left(\psi_{k1k_{2}}^{(i-1)},\cdots,\psi_{kKk_{2}}^{(i-1)}\right)$
$\displaystyle\sim\mathrm{Dir}\left(\epsilon_{0},\cdots,\epsilon_{0}\right),$
$\displaystyle\gamma_{k}^{(i)},c_{k}^{(i)}$
$\displaystyle\sim\mathrm{Gam}\left(\epsilon_{0},\epsilon_{0}\right),\,\mathrm{respectively}.$

The diagrams of three Dirichlet Markov constructions are shown in Figure 3.

##  4 Markov Chain Monte Carlo Inference

In this section, we present a fully-conjugate and efficient Gibbs sampler for
the proposed NS-PGDS. The sampling algorithm depends on two data augmentation
techniques, which we will repeatedly exploit, thus we list them below. Here we
only illustrate the key points of the derivation and the details can be found
in the supplementary material.

######  Lemma 1.

If $y\sim\mathrm{NB}\left(a,g\left(\zeta\right)\right)$ and
$l\sim\mathrm{CRT}\left(y,a\right)$, where $\mathrm{NB}\left(\cdot\right)$
refers to negative binomial distribution, $\mathrm{CRT}\left(\cdot\right)$
represents Chinese restaurant table distribution Teh et al. (2006), and
$g\left(z\right)=1-\mathrm{exp}\left(-z\right)$. Then the joint distribution
of $y$ and $l$ can be equivalently distributed as
$y\sim\mathrm{SumLog}\left(l,g\left(\zeta\right)\right)$ and
$l\sim\mathrm{Pois}\left(a\zeta\right)$ Zhou and Carin (2015), i.e.

$\displaystyle\mathrm{NB}\left(y;a,g\left(\zeta\right)\right)\mathrm{CRT}\left(l;y,a\right)$
$\displaystyle=\mathrm{SumLog}\left(y;l,g\left(\zeta\right)\right)\mathrm{Pois}\left(l;a\zeta\right),$

where $\mathrm{SumLog}\left(l,g\left(\zeta\right)\right)=\sum_{i=1}^{l}x_{i}$
and $x_{i}\sim\mathrm{Log}\left(g\left(\zeta\right)\right)$ are independently
and identically logarithmic distributed random variables Johnson et al.
(2005).

######  Lemma 2.

Suppose $\mathbf{n}=\left(n_{1},\cdots,n_{K}\right)$ and

$\mathbf{n}\mid n\sim\mathrm{DirMult}\left(n,r_{1},\cdots,r_{K}\right),$

where $\mathrm{DirMult}\left(\cdot\right)$ refers to Dirichlet-multimonial
distribution. We sample augmented variable $q\mid
n\sim\mathrm{Beta}\left(n,r_{\cdot}\right)$, where
$r_{\cdot}=\sum_{k=1}^{K}r_{k}$. According to Zhou (2018), conditioning on
$q$, we have $n_{k}\sim\mathrm{NB}\left(r_{k},q\right)$.

Sampling $y_{vk}^{(t)}$ : Use the relationship between Poisson and multinomial
distributions, we sample

$\displaystyle\left(\left(y_{vk}^{(t)}\right)_{k=1}^{K}\mid-\right)\sim\mathrm{Mult}\left(y_{v}^{(t)},\left(\frac{\phi_{vk}\theta_{k}^{(t)}}{\sum_{k=1}^{K}\phi_{vk}\theta_{k}^{(t)}}\right)_{k=1}^{K}\right).$

Sampling $\bm{\phi_{k}}$ : Via Dirichlet-multinomial conjugacy, the posterior
of $\bm{\phi_{k}}$ is

$\displaystyle\left(\bm{\phi_{k}}\mid-\right)\sim\mathrm{Dir}\left(\epsilon_{0}+\textstyle\sum_{t=1}^{T}y_{1k}^{(t)},\cdots,\epsilon_{0}+\textstyle\sum_{t=1}^{T}y_{Vk}^{(t)}\right).$

Sampling $\theta_{k}^{(t)}$ : To sample from the posterior of
$\theta_{k}^{(t)}$, we first sample the auxiliary variables. Set $l_{\cdot
k}^{(T+1)}=0$ and $\zeta^{(T+1)}=0$, sampling backwards from $t=T,\cdots,2$,

$\displaystyle\left(l_{k\cdot}^{(t)}\mid-\right)\sim\mathrm{CRT}\left(y_{\cdot
k}^{(t)}+l_{\cdot
k}^{(t+1)},\tau_{0}\textstyle\sum_{k_{2}=1}^{K}\pi_{kk_{2}}^{i\left(t-1\right)}\theta_{k_{2}}^{(t-1)}\right),$

$\displaystyle\begin{aligned}
&\left(l_{k1}^{(t)},\cdots,l_{kK}^{(t)}\mid-\right)\sim\\\
&\mathrm{Mult}\left(l_{k\cdot}^{(t)},\left(\frac{\pi_{k1}^{i\left(t-1\right)}\theta_{1}^{(t-1)}}{\sum_{k_{2}=1}^{K}\pi_{kk_{2}}^{i\left(t-1\right)}\theta_{k_{2}}^{(t-1)}},\cdots,\frac{\pi_{kK}^{i\left(t-1\right)}\theta_{K}^{(t-1)}}{\sum_{k_{2}=1}^{K}\pi_{kk_{2}}^{i\left(t-1\right)}\theta_{k_{2}}^{(t-1)}}\right)\right).\end{aligned}$

Let us define

$\displaystyle l_{\cdot
k}^{(t)}=\sum_{k_{1}=1}^{K}l_{k_{1}k}^{(t)}\;\mathrm{and}\;\zeta^{(t)}=\mathrm{ln}\left(1+\frac{\delta^{(t)}}{\tau_{0}}+\zeta^{(t+1)}\right).$

After sampling the auxiliary variables, then for $t=1,\cdots,T$, by Poisson-
gamma conjugacy, we obtain

$\displaystyle\left(\theta_{k}^{(1)}\mid-\right)\sim\mathrm{Gam}\left(y_{\cdot
k}^{(1)}+l_{\cdot
k}^{(2)}+\tau_{0}\nu_{k},\tau_{0}+\delta^{(1)}+\zeta^{(2)}\tau_{0}\right),$

$\displaystyle\begin{aligned} &\left(\theta_{k}^{(t)}\mid-\right)\sim\\\
&\mathrm{Gam}\left(y_{\cdot k}^{(t)}+l_{\cdot
k}^{(t+1)}+\tau_{0}\sum_{k_{2}=1}^{K}\pi_{kk_{2}}^{i\left(t-1\right)}\theta_{k_{2}}^{(t-1)},\tau_{0}+\delta^{(t)}+\zeta^{(t+1)}\tau_{0}\right).\end{aligned}$

Sampling $\bm{\Pi}^{(i)}$ : Here we only illustrate Gibbs sampling algorithm
for PR-Gam-Dir construction, sampling algorithms for other constructions can
be found in the supplementary material. We define $M$ as the length of each
interval, and $I$ as the number of intervals. For $i=I,\cdots,2$, because
$\left(l_{1k}^{(i)},\cdots,l_{Kk}^{(i)}\right)$ and $\left(g_{\cdot
1k}^{(i+1)},\cdots,g_{\cdot Kk}^{(i+1)}\right)$ are multinomial distributed,
where $l_{k_{1}k}^{(i)}=\sum_{\left(i-1\right)M+1}^{iM}l_{k_{1}k}^{(t)}$
refers to the summation of $l_{k_{1}k}^{(t)}$ over $i$-th interval. By the
definition of Dirichlet-multinomial distribution and via Lemma 2, defining
$g_{k_{1}k}^{(I+1)}=0$, we sample the auxiliary variables as
$\left(q_{k}^{(i)}\mid-\right)\sim\mathrm{Beta}\left(l_{\cdot
k}^{(i)}+g_{\cdot k}^{(i+1)},\alpha_{\cdot k}^{(i)}\right)$. Via Lemma 2, we
have $\left(l_{k_{1}k}^{(i)}+g_{\cdot
k_{1}k}^{(i+1)}\right)\sim\mathrm{NB}\left(\alpha_{k_{1}k}^{(i)},q_{k}^{(i)}\right)$.
Then we further sample
$\left(h_{k_{1}k}^{(i)}\mid-\right)\sim\mathrm{CRT}\left(l_{k_{1}k}^{(i)}+g_{\cdot
k_{1}k}^{(i+1)},\alpha_{k_{1}k}^{(i)}\right)$. Via Lemma 1, we obtain

$h_{k_{1}k}^{(i)}\sim\mathrm{Pois}\left(-\alpha_{k_{1}k}^{(i)}\mathrm{ln}\left(1-q_{k}^{(i)}\right)\right).$
(5)

For Dirichlet-Randomized-Gamma-Dirichlet Markov construction defined by
Eq.(4), we can equivalently represent it as

$\displaystyle\alpha_{k_{1}k}^{(i)}$
$\displaystyle\sim\mathrm{Gam}\left(g_{k_{1}k}^{(i)}+\epsilon^{\alpha},c_{k}^{(i)}\right),$
(6a) $\displaystyle g_{k_{1}k}^{(i)}$
$\displaystyle=\mathrm{Pois}\left(\gamma^{(i-1)}\textstyle\sum_{k2=1}^{K}\psi^{(i-1)}_{kk_{1}k_{2}}\pi_{k_{2}k}^{(i-1)}\right).$
(6b)

We define
$\lambda_{k_{1}k}^{(i-1)}\triangleq\gamma_{k}^{(i-1)}\sum_{k2=1}^{K}\psi^{(i-1)}_{kk_{1}k_{2}}\pi_{k_{2}k}^{(i-1)}$
for notation conciseness. By Poisson-gamma conjugacy and Eq.(5), we have

$\displaystyle\left(\alpha_{k_{1}k}^{(i)}\mid-\right)\sim\mathrm{Gam}\left(g_{k_{1}k}^{(i)}+\epsilon^{\alpha}+h_{k_{1}k}^{(i)},c_{k}^{(i)}-\mathrm{ln}\left(1-q_{k}^{(i)}\right)\right).$

If $\epsilon^{\alpha}\textgreater 0$, we can sample the posterior of
$g_{k_{1}k}^{(i)}$ via

$\left(g_{k_{1}k}^{(i)}\mid-\right)\sim\mathrm{Bessel}\left(\epsilon^{\alpha}-1,2\sqrt{\alpha_{k_{1}k}^{(i)}c_{k}^{(i)}\lambda_{k_{1}k}^{(i-1)}}\right),$

where $\mathrm{Bessel}\left(\cdot\right)$ denotes Bessel distribution. If
$\epsilon^{\alpha}=0$, we sample $g_{k_{1}k}^{(i)}$ via

$\displaystyle\left(g_{k_{1}k}^{(i)}\mid-\right)\sim\left\\{\begin{array}[]{cc}\mathrm{Pois}\left(\frac{c_{k}^{(i)}\lambda_{k_{1}k}^{(i-1)}}{c_{k}^{(i)}-\mathrm{ln}\left(1-q_{k}^{(i)}\right)}\right)&\mathrm{if}\,h_{k_{1}k}^{(i)}=0\\\
\mathrm{SCH}\left(h_{k_{1}k}^{(i)},\frac{c_{k}^{(i)}\lambda_{k_{1}k}^{(i-1)}}{c_{k}^{(i)}-\mathrm{ln}\left(1-q_{k}^{(i)}\right)}\right)&\mathrm{otherwise},\end{array}\right.$

where $\mathrm{SCH}\left(\cdot\right)$ denotes the shifted confluent
hypergeometric distribution Schein et al. (2019).

Defining $g_{k_{1}k}^{(i)}=g_{k_{1}\cdot
k}^{(i)}=\sum_{k2=1}^{K}g_{k_{1}k_{2}k}^{(i)}$, we first augment

$\displaystyle\left(g_{k_{1}1k}^{(i)},\cdots,g_{k_{1}Kk}^{(i)}\right)\sim\mathrm{Mult}\left(g_{k_{1}k}^{(i)},\left(\psi^{(i-1)}_{kk_{1}k_{2}}\pi_{k_{2}k}^{(i-1)}\right)_{k_{2}=1}^{K}\right).$

Then we obtain
$g_{k_{1}k_{2}k}^{(i)}\sim\mathrm{Pois}\left(\gamma^{(i-1)}\psi^{(i-1)}_{kk_{1}k_{2}}\pi_{k_{2}k}^{(i-1)}\right)$.
By Dirichlet-multinomial conjugacy, we have

$\leavevmode\resizebox{411.93767pt}{}{$\displaystyle\Big{(}\Big{(}\psi_{k1k_{2}}^{(i-1)},\cdots,\psi_{kKk_{2}}^{(i-1)}\Big{)}\mid-\Big{)}\sim\mathrm{Dir}\left(\epsilon_{0}+g_{1k_{2}k}^{(i)},\cdots,\epsilon_{0}+g_{Kk_{2}k}^{(i)}\right)$},$

$\displaystyle\begin{aligned} &\left(\bm{\pi}_{k}^{(i-1)}\mid-\right)\sim\\\
&\mathrm{Dir}\left(\alpha_{1k}^{(i-1)}+l_{1k}^{(i-1)}+g_{\cdot
1k}^{(i)},\cdots,\alpha_{Kk}^{(i-1)}+l_{Kk}^{(i-1)}+g_{\cdot
Kk}^{(i)}\right).\end{aligned}$

Specifically, we have $\alpha_{k_{1}k}^{(1)}=\nu_{k_{1}}\nu_{k}$, if
$k_{1}\neq k$, and $\alpha_{k_{1}k}^{(1)}=\xi\nu_{k}$, if $k_{1}=k$. Via
Poisson-gamma conjugacy, we obtain

$\displaystyle\left(\gamma_{k}^{(i-1)}\mid-\right)\sim\mathrm{Gam}\left(\epsilon_{0}+g_{\cdot
k}^{(i)},\epsilon_{0}+1\right).$

By gamma-gamma conjugacy, we have

$\left(c_{k}^{(i)}\mid-\right)\sim\mathrm{Gam}\left(\epsilon_{0}+\gamma_{k}^{(i-1)},\epsilon_{0}+\textstyle\sum_{k_{1}=1}^{K}\alpha_{k_{1}k}^{(i)}\right).$

##  5 Related Work

Modeling count time sequences has been receiving increasing attentions in
statistical and machine learning communities. Han et al. Han et al. (2014)
adopted linear dynamical systems to capture the underlying dynamics of the
data and leveraged Extended Rank likelihood function to model count
observations. Some Poisson-gamma models assume that the count vector at each
time step is modeled by Poisson factor analysis (PFA) Zhou and Carin (2015)
and leverage special stochastic processes to model the temporal dependencies
of latent factors. For example, gamma process dynamic Poisson factor analysis
(GP-DPFA) Acharya et al. (2015) adopts gamma Markov chains which assumes the
latent factor of the next time step is drawn from a gamma distribution with
the shape parameter be the latent factor of the current time step. Schein et
al. Schein et al. (2016b) proposed Poisson-gamma dynamical systems (PGDSs),
which take the interactions among latent dimensions into account and use a
transition matrix to capture the interactions. Deep dynamic Poisson factor
analysis (DDPFA) Gong and Huang (2017) adopts recurrent neural networks (RNNs)
to capture the complex long-term dependencies of latent factors. Yang and
Koeppl Yang and Koeppl (2018a) applied Poisson-gamma count model to analyze
relational data arising from longitudinal networks, which can capture the
evolution of individual node-group memberships over time. Many modifications
of PGDS have been proposed in recent years. Guo et al. Guo et al. (2018)
proposed deep Poisson-gamma dynamical systems which aim to capture the long-
range temporal dependencies. Schein et al. Schein et al. (2019) employed
Poisson-randomized gamma distribution to build a new transition process of
latent factors. Chen et al. Chen et al. (2021) proposed Switching Poisson-
gamma dynamical systems (SPGDS), allowing PGDS to select from several
transition matrices, and thus can better adapt to nonlinear dynamics. In
contrast to SPGDS, the number of transition matrices of the proposed NS-PGDS
is not limited and thus can be adopted to analyze various complicated non-
stationary count sequences. Filstroff et al. Filstroff et al. (2021)
extensively analyzed many gamma Markov chains for non-negative matrix
factorization and introduced new gamma Markov chains with well-defined
stationary distribution (BGAR).

##  6 Experiments

| GP-DPFA | GMC-RATE | GMC-HIER | BGAR | PGDS |  NS-PGDS (Dir-Dir) |  NS-PGDS (Dir-Gam-Dir) |  NS-PGDS (PR-Gam-Dir)
---|---|---|---|---|---|---|---|---
ICEWS |  MAE | S | 0.259 $\pm 0.005$ | 0.258 $\pm 0.005$ | 0.256 $\pm 0.006$ | 0.264 $\pm 0.006$ | 0.215 $\pm 0.007$ | 0.215 $\pm 0.008$ |  0.214 $\pm 0.008$ | 0.215 $\pm 0.008$
F | 0.176 $\pm 0.005$ | 0.187 $\pm 0.003$ | 0.185 $\pm 0.016$ | 0.222 $\pm 0.043$ | 0.185 $\pm 0.003$ |  0.167 $\pm 0.009$ | 0.169 $\pm 0.006$ | 0.169 $\pm 0.009$
MRE | S | 0.125 $\pm 0.003$ | 0.124 $\pm 0.002$ | 0.122 $\pm 0.003$ | 0.130 $\pm 0.004$ | 0.102 $\pm 0.005$ |  0.101 $\pm 0.005$ |  0.101 $\pm 0.005$ | 0.102 $\pm 0.005$
F | 0.099 $\pm 0.006$ | 0.114 $\pm 0.003$ | 0.111 $\pm 0.018$ | 0.142 $\pm 0.036$ | 0.108 $\pm 0.001$ |  0.094 $\pm 0.005$ | 0.097 $\pm 0.004$ | 0.097 $\pm 0.008$
NIPS |  MAE | S | 18.299 $\pm 6.545$ | 17.105 $\pm 6.449$ | 17.098 $\pm 6.441$ | 17.935 $\pm 6.450$ | 14.706 $\pm 4.414$ | 14.032 $\pm 4.401$ | 14.026 $\pm 4.405$ |  14.014 $\pm 4.387$
F | 48.355 $\pm 1.461$ | 46.234 $\pm 1.629$ | 102.506 $\pm 39.932$ | 62.449 $\pm 14.463$ | 51.562 $\pm 0.679$ |  45.979 $\pm 1.342$ | 46.710 $\pm 1.152$ | 46.582 $\pm 1.196$
MRE | S | 0.729 $\pm 0.412$ | 0.684 $\pm 0.316$ | 0.664 $\pm 0.315$ | 0.769 $\pm 0.366$ | 0.590 $\pm 0.097$ | 0.581 $\pm 0.090$ | 0.581 $\pm 0.090$ |  0.580 $\pm 0.090$
F | 0.415 $\pm 0.016$ |  0.387 $\pm 0.023$ | 0.580 $\pm 0.148$ | 0.465 $\pm 0.049$ | 0.459 $\pm 0.006$ | 0.399 $\pm 0.003$ | 0.395 $\pm 0.006$ | 0.397 $\pm 0.003$
USEI |  MAE | S | 4.681 $\pm 0.564$ | 4.931 $\pm 0.872$ | 4.748 $\pm 0.829$ | 5.244 $\pm 0.939$ | 4.703 $\pm 0.538$ | 4.600 $\pm 0.542$ | 4.608 $\pm 0.541$ |  4.596 $\pm 0.562$
F | 11.665 $\pm 0.367$ | 9.454 $\pm 0.809$ | 12.423 $\pm 1.060$ | 21.948 $\pm 0.133$ | 11.118 $\pm 0.220$ | 7.973 $\pm 1.222$ |  7.168 $\pm 1.221$ | 7.296 $\pm 1.127$
MRE | S | 1.458 $\pm 0.177$ | 1.128 $\pm 0.189$ |  1.088 $\pm 0.162$ | 1.941 $\pm 0.209$ | 1.279 $\pm 0.257$ | 1.309 $\pm 0.220$ | 1.298 $\pm 0.236$ | 1.301 $\pm 0.229$
F | 7.473 $\pm 0.623$ | 6.508 $\pm 0.571$ | 8.929 $\pm 2.514$ | 13.706 $\pm 1.268$ | 4.238 $\pm 0.325$ | 2.602 $\pm 0.455$ |  2.577 $\pm 0.331$ | 2.685 $\pm 0.366$
COVID-19 |  MAE | S | 7.935 $\pm 0.751$ | 7.144 $\pm 1.159$ | 7.240 $\pm 0.848$ | 7.819 $\pm 1.348$ | 7.566 $\pm 1.095$ |  6.969 $\pm 1.107$ | 6.988 $\pm 1.056$ | 6.981 $\pm 1.022$
F | 9.137 $\pm 1.102$ | 9.600 $\pm 1.257$ | 10.409 $\pm 1.910$ | 12.550 $\pm 2.156$ | 9.314 $\pm 0.236$ | 8.799 $\pm 0.706$ |  8.770 $\pm 0.438$ | 9.033 $\pm 0.477$
MRE | S | 0.564 $\pm 0.126$ |  0.493 $\pm 0.136$ | 0.504 $\pm 0.109$ | 0.769 $\pm 0.169$ | 0.558 $\pm 0.130$ | 0.523 $\pm 0.125$ | 0.525 $\pm 0.124$ | 0.526 $\pm 0.123$
F | 0.627 $\pm 0.106$ | 0.556 $\pm 0.052$ | 0.585 $\pm 0.067$ | 0.759 $\pm 0.150$ | 0.585 $\pm 0.007$ | 0.523 $\pm 0.028$ | 0.519 $\pm 0.017$ |  0.513 $\pm 0.014$

Table 1: Results of predictive analysis. ”S” means data smoothing and ”F”
means data forecasting.

We conducted experiments for both predictive and exploratory analysis to
demonstrate the ability of the proposed model in capturing non-stationary
count time sequences. The baseline models included in the experiments are:
$\left.1\right)$ Gamma process dynamic Poisson factor analysis (GP-DPFA)
Acharya et al. (2015). GP-DPFA models the evolution of latent components as
$\theta_{k}^{(t)}\sim\mathrm{Gam}\left(\theta_{k}^{(t-1)},c_{t}\right)$, in
which each component evolves independently of the other components.
$\left.2\right)$ Gamma Markov chains on the rate parameter of gamma
distribution (GMC-RATE) Filstroff et al. (2021). GMC-RATE adopts gamma Markov
chains defined via the rate parameter of the gamma distribution to model the
evolution of $\theta_{k}^{(t)}$ as
$\theta_{k}^{(t)}\sim\mathrm{Gam}\left(\alpha,\beta/\theta_{k}^{(t-1)}\right)$.
$\left.3\right)$ Gamma Markov chains on the rate parameter with hierarchical
auxiliary variable (GMC-HIER) Filstroff et al. (2021). GMC-HIER models the
evolution of latent components with an auxiliary variable as
$z_{k}^{(t)}\sim\mathrm{Gam}\left(\alpha_{z},\beta_{z}\theta_{k}^{(t-1)}\right)$
and
$\theta_{k}^{(t)}\sim\mathrm{Gam}\left(a_{\theta},\beta_{\theta}z_{k}^{(t)}\right)$.
$\left.4\right)$ Autogressive beta-gamma procecss (BGAR) Lewis et al. (1989);
Filstroff et al. (2021). BGAR is also a gamma Markov model. In contrast to the
above models, there is a well-defined stationary distribution for BGAR.
$\left.5\right)$ Poisson-gamma dynamical system (PGDS) Schein et al. (2016b)
takes interactions among latent dimensions into account, and models the
evolution of $\theta_{k}^{(t)}$ as
$\theta_{k}^{(t)}\sim\mathrm{Gam}\left(\tau_{0}\sum_{k_{2}=1}^{K}\pi_{kk_{2}}\theta_{k_{2}}^{(t-1)},\tau_{0}\right)$.

The real-world datasets used in the experiments are: $\left.1\right)$
Integrated Crisis Early Warning System (ICEWS): ICEWS is an international
relations event dataset, comprising interaction events between countries
extracted from news corpora. For ICEWS dataset, we have $T=365$ time steps and
$V=6197$ dimensions, and we set $M=30$. $\left.2\right)$ NIPS: NIPS dataset
contains the papers published in the NeurIPS conference from 1987 to 2015. We
have $T=28$ time steps and $V=2000$ dimensions for NIPS dataset and we set
$M=5$. $\left.3\right)$ U.S. Earthquake Intensity (USEI): USEI dataset
contains a collection of damage and felt reports for U.S. (and a few other
countries) earthquakes. We use the monthly reports from 1957-1986 and have
$T=348,V=64$ and set $M=34$. $\left.4\right)$ COVID-19: COVID-19 dataset
contains daily death cases data for states in the United States, spanning from
March 2020 to June 2020. For this dataset, we have $V=51$ dimensions and
$T=90$ time steps and set $M=20$.

###  6.1 Predictive Analysis

To compare the predictive performance of the proposed model with the
baselines, we considered two standard tasks: data smoothing and forecasting.
For data smoothing task, our objective is to predict $\bm{y}^{(t)}$ given the
remaining data observation $\bm{Y}\backslash\bm{y}^{(t)}$. To this end, we
randomly masked 10 percents of the observed data over non-adjacent time steps,
and predicted the masked values. For forecasting task, we held out data of the
last $S$ time steps, and predicted $\bm{y}^{(T+1)},\cdots,\bm{y}^{(T+S)}$
given $\bm{y}^{(1)},\cdots,\bm{y}^{(T)}$. In this experiment we set $S=2$. We
ran the baseline models including GP-DPFA, PGDS, GMC-RATE, GMC-HIER, BGAR,
using their default settings as provided in Acharya et al. (2015); Schein et
al. (2016b); Filstroff et al. (2021). For NS-PGDS, we set $K=100$ for ICEWS,
$K=10$ for other datasets, and set
$\tau_{0}=1,\gamma_{0}=50,\epsilon_{0}=0.1$. We performed 4000 Gibbs sampling
iterations. In the experiments, we found that the Gibbs sampler started to
converge after 1000 iterations, and thus we set the burn-in time be 2000
iterations. We retained every hundredth sample, and averaged the predictions
over the samples. Mean relative error (MRE) and mean absolute error (MAE) are
adopted to evaluate the model’s predictive capability, which are defined as
$\mathrm{MRE}=\frac{1}{TV}\sum_{t}\sum_{v}\frac{\mid
y_{v}^{(t)}-\hat{y}_{v}^{(t)}\mid}{1+y_{v}^{(t)}}$ and
$\mathrm{MAE}=\frac{1}{TV}\sum_{t}\sum_{v}\mid
y_{v}^{(t)}-\hat{y}_{v}^{(t)}\mid$ respectively, where $y_{v}^{(t)}$ indicates
the true count and $\hat{y}_{v}^{(t)}$ is the prediction.

As the experiment results shown in Table 1, NS-PGDS exhibits improved
performance in both data smoothing and forecasting tasks. We attribute this
enhanced capability to the time-varying transition kernels, which effectively
adapt to the non-stationary environment, and thus achieve improved predictive
performance.

###  6.2 Exploratory Analysis

We used ICEWS and NIPS datasets for exploratory analysis, and chose NS-PGDS
with Dirichlet-Dirichlet Markov chains for illustration. Figure 4(a) and
Figure 4(b) demonstrate the top 2 latent factors inferred by NS-PGDS from
ICEWS dataset. From Figure 4(a) we can see that the main labels are “Iraq
(IRQ)–United States (USA)”, “Iraq (IRQ)–United Kingdom (UK)”, “Russia
(RUS)–United States (USA)”, and so on. This latent factor probably corresponds
to the topic about Iraq war. Besides, in Figure 4(a), there is a peak around
March, 2003, and we know that the Iraq war broke out exactly on 20 March,
2003. In addition, the most dominant labels shown in Figure 4(b) are “Japan
(JPN)–United States (USA)”, “China (CHN)–United States (USA)”, “North Korea
(PRK)–United States (USA)”, “South Korea (KOR)–United States (USA)”, and so
on. We can infer that this latent factor corresponds to “Six-Party Talks” and
other accidents about it.

Figure 4: Latent factors inferred by NS-PGDS. (a) and (b) illustrate the top 2
latent factors inferred from ICEWS dataset, (a) corresponds to Iraq war and
(b) corresponds to the Six-Party Talks. (c) illustrates the evolving trends of
the top 5 latent factors inferred from NIPS dataset.

Figure 4(c) demonstrates the evolving trends of the top 5 latent factors
inferred by NS-PGDS from NIPS dataset, and the legend indicates the
representative words of the corresponding latent factors. Clearly, the green
and blue lines correspond to the latent factors of neural network research
which started to decline from the 1990s. From the 1990s we see that the latent
factors about statistical and probabilistic methods began to dominate the
NeurIPS conference. In addition, NS-PGDS also captured the revival of neural
networks (blue line) from the 2010s. The above observations from the latent
structure inferred by NS-PGDS match our prior knowledge.

Next, we explored the time-varying transition matrices inferred by NS-PGDS. We
chose NIPS dataset for illustratiuon, and set $K=10$ and the interval length
$M$ to be 5. The time-varying transition matrices are shown from Figure 5(b)
to Figure 5(f). At the beginning, matrices shown in Figure 5(b) and Figure
5(c) are close to identity matrices. Then the transition matrices tend to
become block diagonal matrices with 2 blocks, as shown in Figure 5(d)-5(f).
The representative words for latent factors in the first block are “state-
linear-classification”, “network-neural-networks”, “kernel-image-space”,
“network-neural-networks”, “neural-networks-state”. The representative words
for latent factors in the second block are “image-sparse-matrix”, “kernel-
supervised-random”, “matrix-sample-random”, “inference-prior-latent”, “state-
policy-gamma”. The first block primarily captured the correlations among the
research topics about neural networks. The second block reflects that, from
the 1990s, statistical learning and Bayesian methods began to dominate, and
these topics are highly correlated. Figure 5(a) illustrates the transition
matrix inferred by the PGDS, which is averaged over all time steps. Compared
with the NS-PGDS, the PGDS can not capture the informative time-varying
transition dynamics.

Figure 5: Transition matrices inferred from NIPS dataset. (a) illustrates the
transition matrix inferred by the PGDS. (b)-(f) illustrate the time-varying
transition matrices inferred by the NS-PGDS.

We also analyzed the features of the proposed Dirichlet Markov chains. The
left column of Figure 6 demonstrates transition matrices of the first four
sub-intervals of ICEWS dataset inferred by NS-PGDS (Dir-Dir). Because of the
Dir-Dir construction, the consecutive transition matrices smoothly change over
time and thus NS-PGDS may lack sufficient flexibility to capture rapid
dynamics. The middle column of Figure 6 illustrates the transition matrices
inferred by NS-PGDS (Dir-Gam-Dir), which takes mutations among latent
components into account and captured more complicated patterns. Transition
matrices inferred by PR-Gam-Dir construction are shown in the right column of
Figure 6, these matrices not only exhibited sufficient flexibility but also
captured sparser patterns compared with Dir-Gam-Dir construction.

Figure 6: From top to bottom are the first four transition matrices inferred
by different Dirichlet Markov chains from ICEWS dataset. Left column: Matrices
inferred by Dir-Dir construction. Middle column: Matrices inferred by Dir-Gam-
Dir construction. Right column: Matrices inferred by PR-Gam-Dir construction.

##  7 Conclusion

The Poisson-gamma dynamical systems with time-varying transition matrices,
have been proposed to capture complicated dynamics observed in _non-
stationary_ count sequences. In particular, Dirichlet Markov chains are
constructed to allow the underlying transition matrices to evolve over time.
Although the Dirichlet Markov processes lack conjugacy, we have developed
tractable-but-efficient Gibbs sampling algorithms to perform posterior
simulation. The experiment results demonstrate the improved performance of the
proposed NS-PGDS in data smoothing and forecasting tasks, compared with the
PGDS with a stationary transition kernel. Moreover, the experimental results
on several real-world data sets show the explainable structures inferred by
the proposed NS-PGDS. In the future research, we consider to generalize
Dirichlet belief networks by incorporating the proposed Dirichlet Markov chain
constructions, which allow the hierarchical topics to mutate across layers,
and thus can generate more richful text information. We also consider to
capture non-stationary interaction dynamics among individuals over online
social networks Yang and Koeppl (2018c, b, 2020); Yang and Zha (2023, 2024) in
the future research.

## References

  * Acharya et al. [2015] Ayan Acharya, Joydeep Ghosh, and Mingyuan Zhou.  Nonparametric bayesian factor analysis for dynamic count matrices.  In Artificial Intelligence and Statistics, pages 1–9, 2015.
  * Blei and Lafferty [2006] David M Blei and John D Lafferty.  Dynamic topic models.  In Proceedings of the 23rd international conference on Machine learning, pages 113–120, 2006.
  * Chang and Blei [2009] Jonathan Chang and David Blei.  Relational topic models for document networks.  In Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics, pages 81–88, 2009.
  * Chen et al. [2021] Wenchao Chen, Bo Chen, Yicheng Liu, Qianru Zhao, and Mingyuan Zhou.  Switching poisson gamma dynamical systems.  In Proceedings of the Twenty-Ninth International Conference on International Joint Conferences on Artificial Intelligence, pages 2029–2036, 2021.
  * Feng et al. [2021] Luzhao Feng, Ting Zhang, Qing Wang, Yiran Xie, Zhibin Peng, Jiandong Zheng, Ying Qin, Muli Zhang, Shengjie Lai, Dayan Wang, et al.  Impact of covid-19 outbreaks and interventions on influenza in china and the united states.  Nature communications, 12(1):3249, 2021.
  * Filstroff et al. [2021] Louis Filstroff, Olivier Gouvert, Cédric Févotte, and Olivier Cappé.  A comparative study of gamma markov chains for temporal non-negative matrix factorization.  IEEE Transactions on Signal Processing, 69:1614–1626, 2021.
  * Ghahramani and Roweis [1998] Zoubin Ghahramani and Sam T Roweis.  Learning nonlinear dynamical systems using an em algorithm.  In Proceedings of the 11th International Conference on Neural Information Processing Systems, pages 431–437, 1998.
  * Gong and Huang [2017] Chengyue Gong and Win-bin Huang.  Deep dynamic poisson factorization model.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, pages 1665–1673, 2017.
  * Grossman et al. [2020] Guy Grossman, Soojong Kim, Jonah M Rexer, and Harsha Thirumurthy.  Political partisanship influences behavioral responses to governors’ recommendations for covid-19 prevention in the united states.  Proceedings of the National Academy of Sciences, 117(39):24144–24153, 2020.
  * Guo et al. [2018] Dandan Guo, Bo Chen, Hao Zhang, and Mingyuan Zhou.  Deep poisson gamma dynamical systems.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 8451–8461, 2018.
  * Han et al. [2014] Shaobo Han, Lin Du, Esther Salazar, and Lawrence Carin.  Dynamic rank factor model for text streams.  In Proceedings of the 27th International Conference on Neural Information Processing Systems-Volume 2, pages 2663–2671, 2014.
  * Jähnichen et al. [2018] Patrick Jähnichen, Florian Wenzel, Marius Kloft, and Stephan Mandt.  Scalable generalized dynamic topic models.  In International Conference on Artificial Intelligence and Statistics, pages 1427–1435, 2018.
  * Johnson et al. [2005] Norman L Johnson, Adrienne W Kemp, and Samuel Kotz.  Univariate discrete distributions, volume 444.  John Wiley & Sons, 2005.
  * Kalantari and Zhou [2020] Rahi Kalantari and Mingyuan Zhou.  Graph gamma process generalized linear dynamical systems.  arXiv preprint arXiv:2007.12852, 2020.
  * Kalman [1960] R. E. Kalman.  A New Approach to Linear Filtering and Prediction Problems.  Journal of Basic Engineering, 82(1):35–45, 1960.
  * Lewis et al. [1989] Peter AW Lewis, Edward McKenzie, and David Kennedy Hugus.  Gamma processes.  Stochastic Models, 5(1):1–30, 1989.
  * Raymer et al. [2013] James Raymer, Arkadiusz Wiśniowski, Jonathan J Forster, Peter WF Smith, and Jakub Bijak.  Integrated modeling of european migration.  Journal of the American Statistical Association, 108(503):801–819, 2013.
  * Schein et al. [2015] Aaron Schein, John Paisley, David M Blei, and Hanna Wallach.  Bayesian poisson tensor factorization for inferring multilateral relations from sparse dyadic event counts.  In Proceedings of the 21th ACM SIGKDD International conference on knowledge discovery and data mining, pages 1045–1054, 2015.
  * Schein et al. [2016a] Aaron Schein, Mingyuan Zhou, David Blei, and Hanna Wallach.  Bayesian poisson tucker decomposition for learning the structure of international relations.  In International Conference on Machine Learning, pages 2810–2819, 2016.
  * Schein et al. [2016b] Aaron Schein, Mingyuan Zhou, and Hanna Wallach.  Poisson-gamma dynamical systems.  In Proceedings of the 30th International Conference on Neural Information Processing Systems, pages 5012–5020, 2016.
  * Schein et al. [2019] Aaron Schein, Scott W Linderman, Mingyuan Zhou, David M Blei, and Hanna Wallach.  Poisson-randomized gamma dynamical systems.  In Proceedings of the 33rd International Conference on Neural Information Processing Systems, pages 782–793, 2019.
  * Sheldon and Dietterich [2011] Daniel Sheldon and Thomas G Dietterich.  Collective graphical models.  In Proceedings of the 24th International Conference on Neural Information Processing Systems, pages 1161–1169, 2011.
  * Team [2021] IHME COVID-19 Forecasting Team.  Modeling covid-19 scenarios for the united states.  Nature medicine, 27(1):94–105, 2021.
  * Teh et al. [2006] Yee Whye Teh, Michael I Jordan, Matthew J Beal, and David M Blei.  Hierarchical dirichlet processes.  Journal of the American Statistical Association, 101(476):1566–1581, 2006.
  * Unwin et al. [2020] H Juliette T Unwin, Swapnil Mishra, Valerie C Bradley, Axel Gandy, Thomas A Mellan, Helen Coupland, Jonathan Ish-Horowicz, Michaela AC Vollmer, Charles Whittaker, Sarah L Filippi, et al.  State-level tracking of covid-19 in the united states.  Nature Communications, 11(1):1–9, 2020.
  * Wang and McCallum [2006] Xuerui Wang and Andrew McCallum.  Topics over time: a non-markov continuous-time model of topical trends.  In Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 424–433, 2006.
  * Wanner [2021] Philippe Wanner.  How well can we estimate immigration trends using google data?  Quality & Quantity, 55(4):1181–1202, 2021.
  * Wilson [2017] Tom Wilson.  Methods for estimating sub-state international migration: The case of australia.  Spatial Demography, 5(3):171–192, 2017.
  * Winkelmann [2008] Rainer Winkelmann.  Econometric Analysis of Count Data.  Springer Publishing Company, Incorporated, 5th edition, 2008.
  * Yang and Koeppl [2018a] Sikun Yang and Heinz Koeppl.  Dependent relational gamma process models for longitudinal networks.  In International Conference on Machine Learning, pages 5551–5560, 2018.
  * Yang and Koeppl [2018b] Sikun Yang and Heinz Koeppl.  Dependent relational gamma process models for longitudinal networks.  In Proceedings of the International Conference on Machine Learning (ICML), pages 5551–5560, 2018.
  * Yang and Koeppl [2018c] Sikun Yang and Heinz Koeppl.  A Poisson gamma probabilistic model for latent node-group memberships in dynamic networks.  In Proceedings of the AAAI Conference on Artificial Intelligence (AAAI), pages 4366–4373, 2018.
  * Yang and Koeppl [2020] Sikun Yang and Heinz Koeppl.  The Hawkes edge partition model for continuous-time event-based temporal networks.  In Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI), pages 460–469, 2020.
  * Yang and Zha [2023] Sikun Yang and Hongyuan Zha.  Estimating latent population flows from aggregated data via inversing multi-marginal optimal transport.  In Proceedings of the 2023 SIAM International Conference on Data Mining (SDM), pages 181–189, 2023.
  * Yang and Zha [2024] Sikun Yang and Hongyuan Zha.  A variational autoencoder for neural temporal point processes with dynamic latent graphs.  In Proceedings of the AAAI Conference on Artificial Intelligence (AAAI), 2024.
  * Yuan and Kalbfleisch [2000] Lin Yuan and John D Kalbfleisch.  On the bessel distribution and related problems.  Annals of the Institute of Statistical Mathematics, 52:438–447, 2000.
  * Zhou and Carin [2012] M Zhou and L Carin.  Augment-and-conquer negative binomial processes.  Advances in Neural Information Processing Systems, 4:2546–2554, 2012.
  * Zhou and Carin [2015] Mingyuan Zhou and Lawrence Carin.  Negative binomial process count and mixture modeling.  IEEE Transactions on Pattern Analysis & Machine Intelligence, 37(02):307–320, 2015.
  * Zhou [2018] Mingyuan Zhou.  Nonparametric bayesian negative binomial factor analysis.  Bayesian Analysis, 13(4):1065–1093, 2018.


:
2022 Conference on Health, Inference, and Learning (CHIL) 2022

# Uncertainty-Aware Text-to-Program for
Question Answering on Structured Electronic Health Records

Daeyoung Kim<EMAIL_ADDRESS>
Seongsu Bae<EMAIL_ADDRESS>
Seungho Kim<EMAIL_ADDRESS>
Edward Choi<EMAIL_ADDRESS>
KAIST  Republic of Korea

###### Abstract

sections/0_abstract

### Data and Code Availability

Our source code and dataset are available on the official
repository111https://github.com/cyc1am3n/text2program-for-ehr.

##  1 Introduction

sections/1_introduction

##  2 Related Works

sections/2_related_works

Figure 1:  Illustrative examples of the natural language question (NLQ) and
the corresponding program, composed of several predefined operations. To
answer a single natural language question, we have to execute a series of
atomic operations in sequence.

tab:grammar  Operation Description Return Data Type gen_entset_down(entSet,
rel) the set of object entities associated with relation rel for each subject
entity in the entSet entSet gen_entset_up(rel, entSet) the set of subject
entities associated with relation rel for each object entity in the entSet
entSet gen_litset(entSet, rel) the set of literal values associated with the
relation rel for each subject entity of entSet litSet gen_entset_equal(rel,
lit) the set of subject entities which have literal value equal to lit for
relation rel entSet gen_entset_atleast(rel, lit) the set of subject entities
which have literal value of at least lit for relation rel entSet
gen_entset_atmost(rel, lit) the set of subject entities which have literal
value of at most lit for relation rel entSet gen_entset_less(rel, lit) the set
of subject entities which have smaller literal value than lit for relation rel
entSet gen_entset_more(rel, lit) the set of subject entities which have
greater literal value than lit for relation rel entSet count_entset(entSet)
the number of entities in entSet int intersect_entsets(entSet1, entSet2) the
set of entities that exist in common in entSet1 and entSet2 entSet
maximum_litset(litSet) the largest value in litSet float
minimum_litset(litSet) the smallest value in litSet float
average_litset(litSet) the average value of litSet float
concat_litsets(litSet1, litSet2) the combined list of litSet1 and litSet2
litSets

Table 1: Description of the custom-defined operations and return data types

##  3 Methodology

sections/3_methodology

##  4 Experiments

sections/4_experiments

##  5 Conclusion

sections/5_conclusion

## Institutional Review Board (IRB)

This research does not require IRB approval.

This work was supported by Institute of Information & Communications
Technology Planning & Evaluation (IITP) grant (No.2019-0-00075, Artificial
Intelligence Graduate School Program(KAIST)) and National Research Foundation
of Korea (NRF) grant (NRF-2020H1D3A2A03100945), funded by the Korea government
(MSIT).

## References

  * Bae et al. (2021) Seongsu Bae, Daeyoung Kim, Jiho Kim, and Edward Choi.  Question answering for complex electronic health records database using unified encoder-decoder architecture.  In _Machine Learning for Health_ , pages 13–25. PMLR, 2021.
  * Chang et al. (2020) Jie Chang, Zhonghao Lan, Changmao Cheng, and Yichen Wei.  Data uncertainty learning in face recognition.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 5710–5719, 2020.
  * Der Kiureghian and Ditlevsen (2009) Armen Der Kiureghian and Ove Ditlevsen.  Aleatory or epistemic? does it matter?  _Structural safety_ , 31(2):105–112, 2009.
  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_ , 2018.
  * Dusenberry et al. (2020) Michael W Dusenberry, Dustin Tran, Edward Choi, Jonas Kemp, Jeremy Nixon, Ghassen Jerfel, Katherine Heller, and Andrew M Dai.  Analyzing the role of model uncertainty for electronic health records.  In _Proceedings of the ACM Conference on Health, Inference, and Learning_ , pages 204–213, 2020.
  * Hua et al. (2020) Yuncheng Hua, Yuan-Fang Li, Guilin Qi, Wei Wu, Jingyao Zhang, and Daiqing Qi.  Less is more: Data-efficient complex question answering over knowledge bases.  _Journal of Web Semantics_ , 65:100612, 2020.
  * Johnson et al. (2016) Alistair EW Johnson, Tom J Pollard, Lu Shen, H Lehman Li-Wei, Mengling Feng, Mohammad Ghassemi, Benjamin Moody, Peter Szolovits, Leo Anthony Celi, and Roger G Mark.  Mimic-iii, a freely accessible critical care database.  _Scientific data_ , 3(1):1–9, 2016.
  * Kendall and Gal (2017) Alex Kendall and Yarin Gal.  What uncertainties do we need in bayesian deep learning for computer vision?  _Advances in neural information processing systems_ , 30, 2017.
  * Lakshminarayanan et al. (2017) Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell.  Simple and scalable predictive uncertainty estimation using deep ensembles.  _Advances in neural information processing systems_ , 30, 2017.
  * Liang et al. (2017) Chen Liang, Jonathan Berant, Quoc Le, Kenneth Forbus, and Ni Lao.  Neural symbolic machines: Learning semantic parsers on freebase with weak supervision.  In _Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 23–33, 2017.
  * Lin (2004) Chin-Yew Lin.  Rouge: A package for automatic evaluation of summaries.  In _Text summarization branches out_ , pages 74–81, 2004.
  * Luong et al. (2015) Minh-Thang Luong, Hieu Pham, and Christopher D Manning.  Effective approaches to attention-based neural machine translation.  In _Proceedings of the 2015 Conference on Empirical Methods in Natural Language Processing_ , pages 1412–1421, 2015.
  * Malinin and Gales (2020) Andrey Malinin and Mark Gales.  Uncertainty estimation in autoregressive structured prediction.  In _International Conference on Learning Representations_ , 2020.
  * Ott et al. (2018) Myle Ott, Michael Auli, David Grangier, and Marc’Aurelio Ranzato.  Analyzing uncertainty in neural machine translation.  In _International Conference on Machine Learning_ , pages 3956–3965. PMLR, 2018.
  * Pampari et al. (2018) Anusri Pampari, Preethi Raghavan, Jennifer Liang, and Jian Peng.  emrqa: A large corpus for question answering on electronic medical records.  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing_ , pages 2357–2368, 2018.
  * Park et al. (2021) Junwoo Park, Youngwoo Cho, Haneol Lee, Jaegul Choo, and Edward Choi.  Knowledge graph-based question answering with electronic health records.  In _Machine Learning for Healthcare Conference_ , pages 36–53. PMLR, 2021.
  * Pollard et al. (2018) Tom J Pollard, Alistair EW Johnson, Jesse D Raffa, Leo A Celi, Roger G Mark, and Omar Badawi.  The eicu collaborative research database, a freely available multi-center database for critical care research.  _Scientific data_ , 5(1):1–13, 2018.
  * Raffel et al. (2020) Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J Liu.  Exploring the limits of transfer learning with a unified text-to-text transformer.  _Journal of Machine Learning Research_ , 21:1–67, 2020\.
  * Raghavan et al. (2021) Preethi Raghavan, Jennifer J Liang, Diwakar Mahajan, Rachita Chandra, and Peter Szolovits.  emrkbqa: A clinical knowledge-base question answering dataset.  In _Proceedings of the 20th Workshop on Biomedical Language Processing_ , pages 64–73, 2021.
  * Saha et al. (2018) Amrita Saha, Vardaan Pahuja, Mitesh M Khapra, Karthik Sankaranarayanan, and Sarath Chandar.  Complex sequential question answering: Towards learning to converse over linked question answer pairs with a knowledge graph.  In _Thirty-Second AAAI Conference on Artificial Intelligence_ , 2018\.
  * Saha et al. (2019) Amrita Saha, Ghulam Ahmed Ansari, Abhishek Laddha, Karthik Sankaranarayanan, and Soumen Chakrabarti.  Complex program induction for querying knowledge bases in the absence of gold programs.  _Transactions of the Association for Computational Linguistics_ , 7:185–200, 2019.
  * Suster and Daelemans (2018) Simon Suster and Walter Daelemans.  Clicr: a dataset of clinical case reports for machine reading comprehension.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ , pages 1551–1563, 2018.
  * Wang et al. (2020) Ping Wang, Tian Shi, and Chandan K Reddy.  Text-to-sql generation for question answering on electronic medical records.  In _Proceedings of The Web Conference 2020_ , pages 350–361, 2020\.
  * Wiseman and Rush (2016) Sam Wiseman and Alexander M Rush.  Sequence-to-sequence learning as beam-search optimization.  In _Proceedings of the 2016 Conference on Empirical Methods in Natural Language Processing_ , pages 1296–1306, 2016.
  * Xiao and Wang (2021) Yijun Xiao and William Yang Wang.  On hallucination and predictive uncertainty in conditional language generation.  In _Proceedings of the 16th Conference of the European Chapter of the Association for Computational Linguistics: Main Volume_ , pages 2734–2744, 2021.
  * Xu et al. (2020) Jiacheng Xu, Shrey Desai, and Greg Durrett.  Understanding neural abstractive summarization models via uncertainty.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 6275–6281, 2020.
  * Yih et al. (2015) Scott Wen-tau Yih, Ming-Wei Chang, Xiaodong He, and Jianfeng Gao.  Semantic parsing via staged query graph generation: Question answering with knowledge base.  In _Proceedings of the Joint Conference of the 53rd Annual Meeting of the ACL and the 7th International Joint Conference on Natural Language Processing of the AFNLP_ , 2015.

sections/6_appendix