\section{Introduction}
\label{S1:Introduction}

Predicting extreme volatility events in financial markets is essential
when estimating risk. A standard approach to extreme event prediction is
to find the precursory patterns prior to an extreme event or to quantify
the probability that a given pattern is a precursor to an extreme event
\citep{Hallerberg-Altmann-Holstein-Kantz-2007-PRE,Hallerberg-Kantz-2008-PRE}.
\cite{Bogachev-Bunde-2009-PRE,Bogachev-Bunde-2011-PA} propose a new method based on the statistics of the recurrence
intervals between events exceeding a threshold to determine the risk
probability $W(\Delta t | t)$ that an extreme event will occur within
the next $\Delta t$ intervals when the last extreme event occurred $t$
periods ago. They
find that when examining real market data and model data with a low
level of noise the predicting method based on recurrence interval
analysis produces a better performance forecast than the method based on
precursor pattern recognition
\citep{Bogachev-Bunde-2009-PRE,Bogachev-Bunde-2011-PA}.

Understanding of the recurrence interval, defined as the waiting time
between consecutive events with values greater than a predefined
threshold $Q$, is essential in uncovering the underlying laws governing
extreme events in many fields. Recurrent interval analysis has been
carried out on many kinds of time series in predicting the probability
that an extreme event will occur, including records of climate
\citep{Bunde-Eichner-Havlin-Kantelhardt-2004-PA,Bunde-Eichner-Kantelhardt-Havlin-2005-PRL},
seismic activities \citep{Saichev-Sornette-2006-PRL}, energy dissipation
rates of three-dimensional turbulence
\citep{Liu-Jiang-Ren-Zhou-2009-PRE}, heartbeat intervals in medical
science \citep{Bogachev-Kireenkov-Nifontov-Bunde-2009-NJP},
precipitation and river runoff \citep{Bogachev-Bunde-2012-EPL}, internet
traffic \citep{Bogachev-Bunde-2009-EPL,Cai-Fu-Zhou-Gu-Zhou-2009-EPL},
financial volatilities
\citep{Yamasaki-Muchnik-Havlin-Bunde-Stanley-2005-PNAS,Xie-Jiang-Zhou-2014-EM},
equity returns
\citep{Yamasaki-Muchnik-Havlin-Bunde-Stanley-2006-inPFE,Bogachev-Eichner-Bunde-2007-PRL,Bogachev-Bunde-2008-PRE,Bogachev-Bunde-2009-PRE,Ren-Zhou-2010-NJP,Ludescher-Tsallis-Bunde-2011-EPL,He-Chen-2011b-PA,Meng-Ren-Gu-Xiong-Zhang-Zhou-Zhang-2012-EPL,Suo-Wang-Li-2005-EM},
and trading volumes
\citep{Podobnik-Horvatic-Petersen-Stanley-2009-PNAS,Ren-Zhou-2010-PRE,Li-Wang-Havlin-Stanley-2011-PRE}.
An improved method of estimating value at risk (VaR) in financial
markets has been proposed based on the recurrence interval between the
last two returns below $-Q$. This method is significantly more accurate
than traditional estimates based on the overall or local return
distributions
\citep{Bogachev-Bunde-2009-PRE,Ludescher-Tsallis-Bunde-2011-EPL}.

To accurately estimate the risk probability and the VaR based on
recurrence interval analysis we need the set of distribution and
memory behavior of the recurrence time between extreme events. It is found that
the recurrence intervals of the time series in many different fields exhibit fat-tailed distributions and long- and short-range memories,
indicating that extreme events are not described by the Poisson
process. Unlike long- and short-range memory behaviors, which are easily
testable using conditional distribution analysis and the DFA method, the
distribution form of recurrence intervals is still elusive.  For
example, in financial markets the recurrence intervals of different data
types (return, volatility, and trading volume), different data
resolutions (minute-by-minute and daily), and different markets fit
different distributions, including power-law, stretched exponential, and
$q$-exponential.  It has been found that the recurrence interval
distribution above a fixed threshold has a power-law tail for the daily
volatilities in Japanese market
\citep{Kaizoji-Kaizoji-2004a-PA,Yamasaki-Muchnik-Havlin-Bunde-Stanley-2005-PNAS},
the minute-by-minute volatilities in Korean market
\citep{Lee-Lee-Rikvold-2006-JKPS} and Italian market
\citep{Greco-SorrisoValvo-Carbone-Cidone-2008-PA}, the daily returns in
US stock markets
\citep{Bogachev-Eichner-Bunde-2007-PRL,Bogachev-Bunde-2008-PRE,Bogachev-Bunde-2009-PRE},
the minute-by-minute returns in Chinese markets
\citep{Ren-Zhou-2010-NJP}, and the minute-by-minute trading volume in US
markets \citep{Li-Wang-Havlin-Stanley-2011-PRE} and Chinese markets
\citep{Ren-Zhou-2010-PRE}. A number of studies ranging from daily to
high-frequency data and from developed to emerging markets
\citep{Wang-Yamasaki-Havlin-Stanley-2006-PRE,Wang-Weber-Yamasaki-Havlin-Stanley-2007-EPJB,Jung-Wang-Havlin-Kaizoji-Moon-Stanley-2008-EPJB,Qiu-Guo-Chen-2008-PA,Ren-Gu-Zhou-2009-PA,Ren-Guo-Zhou-2009-PA,Jeon-Moon-Oh-Yang-Jung-2010-JKPS,Wang-Wang-2012-CIE,Xie-Jiang-Zhou-2014-EM},
have also reported that the distribution of the recurrence intervals of
financial volatility is a stretched
exponential. Reference~\cite{Suo-Wang-Li-2005-EM} reports that in
Chinese markets the recurrence time between returns above a given
positive threshold or below a negative threshold for the index spot and
futures fits a stretched exponential distribution.  The recurrence
intervals between losses in financial markets were recently found to fit
a $q$-exponential distribution
\citep{Ludescher-Tsallis-Bunde-2011-EPL,Ludescher-Bunde-2014-PRE}.

In this paper we describe the datasets in Sec.~\ref{S1:Data}, present
the theoretical framework for predicting large volatilities in
Sec.~\ref{S1:Framework}, determine the distribution of the recurrence
intervals between large volatilities in Sec.~\ref{S1:Distribution}, and
report the hazard probability results and predicting algorithm
performance in Sec.~\ref{S1:Risk}.  In Sec.~\ref{S1:Conclusion} we
summarize our findings.

\section{Data description}
\label{S1:Data}

To carry out a detailed recurrence interval analysis of Chinese stock
markets, we include as many Chinese stocks in our analyzing sample as
possible. The minute-by-minute price data of all stocks in the Chinese
markets are extracted from the RESSET financial database. The extracting
period is from 26 July 1999 to 30 December 2011, which is the maximum
spanning period allowed in the RESSET database. To ensure that the
recurrent interval results between the top 1\% volatilities will have
more than 1,000 data points, we select only those stocks that have a
minimum of two years of trading records. Having this large a sample size
lowers the error rate when we use a maximum likelihood estimation to fit
the distributions.  Finally, we have 1891 stocks in our sample, which include 853 A-shares,
54 B-shares, and 63 ChiNext shares in the Shenzhen market, and 867
A-shares and 54 B-shares in the Shanghai market.

\section{Framework of predicting large volatilities}
\label{S1:Framework}
\subsection{Hazard probability $W(\Delta t | t)$}

We use the hazard probability $W(\Delta t | t)$ to forecast the
occurrance of large volatility events. The $W(\Delta t | t)$ is the
probability that there will be additional waiting time $\Delta t$ before
another large volatility event occurs when the previous large volatility
event occurred $t$ time ago. This probability is the key early-warning
measurement for the occurrence of extreme volatilities. The early
warning is triggered when the probability $W(\Delta t | t)$ is greater
than a predefined alarm threshold. We can theoretically derive this
hazard probability if we have the distribution of the time intervals
between consecutive extreme volatilities, which are defined as the
volatilities that exceed a given threshold $Q$.

Using the probability density $p(t)$ of the recurrence intervals between
the extreme volatilities, if the time elapsed since the last extreme
event is $t$, we want to determine the probability density function
$p(\Delta{t}|t)$ that quantifies the additional waiting time $\Delta{t}$
until the next extreme event. Using the Bayes theorem for conditional
probabilities, the probability that an event A occurs, given the
knowledge of an event B, is simply the quotient of the probability of
the event A without constraint and the probability of event B
\citep{Sornette-Knopoff-1997-BSSA},
\begin{equation}
  p(A|B) = \frac{p(AB)}{p(B)},
  \label{Eq:P:A:B}
\end{equation}
where $p(AB)=p(t+\Delta{t})$ is the probability that no event occurs
from 0 to $t$ and that an event occurs at $t+\Delta{t}$, and
$p(B)=\int_t^{\infty}p(s)ds$ is the probability that no event occurs
from 0 to $t$. Thus we have
\begin{equation}
  p(\Delta{t}|t) = \frac{p(t+\Delta{t})}{\int_t^{\infty}p(s)ds}.
  \label{Eq:P:Dt:t}
\end{equation}
Thus the hazard probability $W(\Delta{t},t)$ that an extreme event will
occur after a short time $\Delta{t}\ll t$ since the occurrance of the
previous extreme event can be expressed
\begin{equation}
 W(\Delta{t}|t)=\frac{\int_t^{t+\Delta{t}}p(t)dt}{\int_t^{\infty}p(t)dt}.
 \label{Eq:Wq}
\end{equation}

\subsection{Predicting algorithm}

For a given distribution $p(t)$ of the recurrence intervals between
extreme volatilities, the formula for $W(\Delta{t}|t)$ can be obtained
using equation~(\ref{Eq:Wq}). If $\Delta t=1$ for $W(\Delta{t}|t)$, the
hazard probability and a decision-making algorithm can be used to
predict large volatilities \citep{Bogachev-Bunde-2011-PA}. To trigger an
early warning that a large volatility is about to occur, we set a
threshold $Q_p$ for the hazard probability. When the hazard probability
exceeds $Q_p$, an alarm that a large volatility will occur during the
next time point is activated. We next estimate the $Q_p$ parameter,
which is the maximum correct prediction rate when the maximum
false-alarm tolerance is set.

To determine $Q_p$ we estimate the correct prediction and false alarm
rates for each $Q_p$ in the range of $[0, 1]$. We then plot the correct
prediction rate with respect to the false alarm rate and get the
``receiver operator characteristic'' (ROC) curve
\citep{Bogachev-Bunde-2009-PRE,Bogachev-Bunde-2009-EPL,Bogachev-Kireenkov-Nifontov-Bunde-2009-NJP,Bogachev-Bunde-2011-PA}.
The ROC curve is used to quantify prediction efficiency. The satisfied
$Q_p$ corresponds to the point at which the false alarm rate equals the
tolerant alarm level on the ROC curve.

To estimate the correct prediction and false alarm rates we generate for
a given $Q_p$ two forecasting signals---alarms and non-alarms---at each
time point.  By comparing the forecasting signals with the real data, we
obtain one of four outcomes at each time point
\citep{Bogachev-Bunde-2011-PA}, (i) a correct prediction of a large
volatility event, (ii) a correct prediction of a non-large volatility
event, (iii) a missed event, and (iv) a false alarm. By recording in our
testing records how many times each outcome occurs we can estimate the
correct prediction rate $D$ and the false alarm rate $A$ using
\begin{equation}
 D = \frac{O_{11}}{ O_{01} + O_{11}},~~A = \frac{O_{10}}{O_{00} + O_{10}},
 \label{Eq:Wq:AD}
\end{equation}
where $O_{11}$ is the number of large volatility events that are
correctly predicted, $O_{00}$ the number of non-large volatility events
that are correctly predicted, $O_{01}$ the number of missed events, and
$O_{10}$ the number of false alarms. All possible pairs of $(D, A)$ will
be obtained if we vary the $Q_p$ range from 0 to 1.

By definition the ROC curve will be $D=A=1$ if $Q_p =0$ and $D=A=0$ if
$Q_p = 1$. Note that the ROC curve joins the point $(0, 0)$ in the left
bottom corner to the point $(1, 1)$ in the right top corner. Note also
that, for the random guess outcome, $D=A$, a straight line between the two
corners. This occurs when there is no memory in the data. For a fixed
value of $A$, the larger the value of the correct prediction rate, the
better this algorithm performs.

\section{Distribution of recurrence intervals between large volatilities}
\label{S1:Distribution}

\subsection{Definition of volatilities and recurrence intervals}

For a given minute-by-minute price series $p(t)$, the minute-by-minute
volatility $\omega(t)$ can be estimated \citep{Bollen-Inder-2002-JEF}
using
\begin{equation}
  \omega(t)=|\ln p(t)-\ln p(t-1)|.
  \label{vlt}
\end{equation}
In order to eliminate the influence of the daily periodic patterns, we
remove the intraday patterns from the volatility series $\omega(t)$ on
each trading day,
\begin{equation}
  \omega'(s)=\omega(s)/A(s),
  \label{vltrip}
\end{equation}
where $A(s) = \sum_i^N \omega(i, s) / N$. Here $\omega(i, s)$ represents
the volatility at time $s$ on day $i$. The normalized volatility series
$v(t)$ is then obtained by dividing $\omega'(t)$ by its standard
deviation,
\begin{equation}
  v(t)= \frac{\omega'(t)}{ \sqrt{[ \langle \omega'(t)^2 \rangle - \langle \omega'(t) \rangle^2 ]} }.
  \label{nvlt}
\end{equation}

The focus of our study is the recurrence interval between the normalized
volatilities exceeding a predefined threshold $Q$. To compare the
results between different stocks, we quantify $Q$ by its mean recurrence
time $\tau_Q$. There is a one-to-one correspondence between $Q$ and
$\tau_Q$, such that \citep{Podobnik-Horvatic-Petersen-Stanley-2009-PNAS}
\begin{equation}
  \frac{1}{\tau_Q} = \int_Q^\infty p(v) {\rm{d}}v,
  \label{tauQpv}
\end{equation}

where $p(v)$ is the probability distribution of the volatility. Here we
restrict $\tau_Q$ to a range of $[20, 100]$. This range corresponds to
the extreme volatilities from a top value of \%5 to 1\%, which is often
considered in the risk estimation.

\subsection{Distribution formula of recurrence intervals}

To analytically determine the hazard probability we find the
distribution that best approximates the recurrence interval distribution
for all the stocks in our sample. We also determine whether the
distribution parameters are dependent on the mean recurrence interval
$\tau_Q$ and on whether the market is bear or bull. Previous research
has indicated that the distribution of recurrence intervals between
returns below a negative threshold $-Q$ depends only on the mean
recurrence interval $\tau_Q$, and not on a specific asset or on the time
resolution of the data
\citep{Ludescher-Tsallis-Bunde-2011-EPL,Ludescher-Bunde-2014-PRE}. We
first check whether the recurrence time between volatilities above a
threshold $Q$ exhibits this behavior. Figure~\ref{Fig:RI:PDF}(a) shows
the probability distribution of the recurrence intervals for 10 randomly
chosen stocks. We transform the distribution curves of different
$\tau_Q$ values by a factor to increase visibility. Note that for the
same $\tau_Q$ value the recurrence time distributions of different
stocks nearly overlap on the same curve. This means that the return
intervals between volatilities that exceed a threshold $Q$ may exhibit a
universal distribution for different stocks when $\tau_Q$ is fixed. We
also want to know whether the distributions of different $\tau_Q$ values
share the same pattern. Previous research indicates that these
distributions are influenced only by the mean recurrence time $\tau_Q$
for return recurrence intervals
\citep{Ludescher-Bunde-2014-PRE}. Although the distributions in
Fig.~\ref{Fig:RI:PDF}(a) seem to be different for different $\tau_Q$, if
we scale the distribution using the mean recurrence time $\tau_Q$ the
six distributions seem to be parallel [see Fig.~\ref{Fig:RI:PDF}(b)].

\begin{figure}[htb]
\centering
\includegraphics[width=8cm]{Fig_PDF_RI_Stocks.eps}
\includegraphics[width=7.3cm]{Fig_PDF_ScaledRI_Stocks.eps}
\caption{\label{Fig:RI:PDF} (color online). Probability distribution of
  the recurrence intervals for 10 randomly chosen stocks. For better
  visibility, the distribution curve of $\tau_Q =25$, $40$, $60$, $80$,
  and $100$ are shifted vertically by a factor of 10, 100, 1000, 10000,
  and 100000, respectively.  (a) Original recurrence intervals. (b)
  Scaled recurrence intervals.}
\end{figure}

To quantitatively measure how the recurrence interval
distributions vary with the mean interval $\tau_Q$, we need a
suitable distributional formula for capturing the recurrence interval
distribution. Previous research shows that the recurrence intervals can be fitted
by the stretched distribution
\citep{Wang-Yamasaki-Havlin-Stanley-2006-PRE,Wang-Weber-Yamasaki-Havlin-Stanley-2007-EPJB,Jung-Wang-Havlin-Kaizoji-Moon-Stanley-2008-EPJB,Qiu-Guo-Chen-2008-PA,Jeon-Moon-Oh-Yang-Jung-2010-JKPS,Wang-Wang-2012-CIE,Suo-Wang-Li-2005-EM},
the power-law distribution with an exponential cutoff
\citep{Kaizoji-Kaizoji-2004a-PA,Yamasaki-Muchnik-Havlin-Bunde-Stanley-2005-PNAS,Lee-Lee-Rikvold-2006-JKPS,Greco-SorrisoValvo-Carbone-Cidone-2008-PA},
and the $q$-exponential distribution
\citep{Ludescher-Tsallis-Bunde-2011-EPL,Ludescher-Bunde-2014-PRE}.
The Weibull distribution is usually comparable to the $q$-exponential distribution
\citep{Poloti-Scalas-2008-PA,Jiang-Chen-Zhou-2008-PA}, and we use both
two-parameter and three-parameter Weibull distributions to fit
the recurrence intervals in our analysis. The following are five
candidate distributions: the stretched exponential distribution,
\begin{equation}
 p(\tau) = a \exp[-(b \tau)^\mu],
 \label{Eq:StrExp:Tau:PDF}
\end{equation}
the power-law distribution with an exponential cutoff,
\begin{equation}
 p(\tau) = c \tau^{-\gamma-1} \exp(-k \tau),
 \label{Eq:PowExp:Tau:PDF}
\end{equation}
the $q$-exponential distribution,
\begin{equation}
 p(\tau) = (2-q) \lambda [1+(q-1) \lambda \tau]^{-\frac{1}{q-1}},
 \label{Eq:qExp:Tau:PDF}
\end{equation}
the two-parameter Weibull distribution,
\begin{equation}
 p(\tau) = \frac{\zeta}{d} \left( \frac{\tau}{d}\right)^{\zeta-1}\exp
 \left[-\left( \frac{\tau}{d}
\right)^\zeta \right],
 \label{Eq:2WBL:Tau:PDF}
\end{equation}
and the three-parameter Weibull distribution,
\begin{equation}
 p(\tau) =\frac{\zeta}{d} \left( \frac{\tau -
   \tau_0}{d}\right)^{\zeta-1}\exp \left[- \left( \frac{\tau -
     \tau_0}{d}\right)^\zeta  \right].
 \label{Eq:3WBL:Tau:PDF}
\end{equation}

To compare the recurrence interval distributions of different $\tau_Q$
values, we normalize the recurrence intervals by $\tau_Q$, i.e., $x =
\tau / \tau_Q$. We obtain the five candidate distributions used to fit
the normalized recurrence intervals by substituting $\tau = x \tau_Q$
into $f(x) = p(\tau) \tau_Q$, i.e.,
\begin{equation}
f(x) = a \bar{\tau} \exp[-(b \tau_Q x)^\mu],
 \label{Eq:StrExp:x:PDF}
\end{equation}
\begin{equation}
f(x) = c \tau_Q^{-\gamma} x^{-\gamma-1} \exp(-k \tau_Q x),
 \label{Eq:PowExp:x:pdf}
\end{equation}
\begin{equation}
 p(\tau) = (2-q) (\lambda \tau_Q) [1+(q-1) (\lambda \tau_Q)
   \tau]^{-\frac{1}{q-1}},
 \label{Eq:qExp:x:PDF}
\end{equation}
\begin{equation}
 f(x) = \frac{\zeta}{d/\tau_Q} \left(
 \frac{x}{d/\tau_Q}\right)^{\zeta-1}\exp\left[-\left(
   \frac{x}{d/\tau_Q}\right)^\zeta \right],
 \label{Eq:2WBL:x:pdf}
\end{equation}
and
\begin{equation}
f(x) =\frac{\zeta}{d/\tau_Q} \left( \frac{x -
  \tau_0/\tau_Q}{d/\tau_Q}\right)^{\zeta-1}\exp\left[-\left( \frac{x -
    \tau_0/\tau_Q}{d/\tau_Q}\right)^\zeta \right].
 \label{Eq:3WBL:x:pdf}
\end{equation}

We use the maximum likelihood estimation (MLE) method to estimate the
parameters of the five distribution parameters. The details of fitting
the stretched exponential distribution and the power-law distribution
with an exponential cutoff are presented in the Appendix. Note that in
the following analysis we fit only the scaled recurrence intervals.

\begin{figure}[htb]
\centering
\includegraphics[width=8cm]{Fig_PDF_ScaledRI_Fit_000001.eps}
\includegraphics[width=8cm]{Fig_PDF_ScaledRI_Fit_900956.eps}
\caption{\label{Fig:RI:PDF:Fit} (color online). Plots of distribution
  fits to normalized recurrence intervals for two stocks. The open
  markers are the empirical distributions and the solid lines are the
  fits to the five candidate distributions. For better visibility, the
  curves of $\tau_Q =25$, $40$, $60$, $80$, and $100$ are shifted
  vertically by a factor of 10, 100, 1000, 10000, and 100000,
  respectively.  (a) Stock 000001. (b) Stock 900956.}
\end{figure}

Figure~\ref{Fig:RI:PDF:Fit} shows the empirical distributions and the
fitting results of the five candidate distributions of the scaled
recurrence intervals for two stocks, 000001 and 900956. Note that in the
central regions of the distributions all five candidates agree with the
empirical data. Note also that the $q$-exponential distribution better
fits the distribution tail for all $\tau_Q$ than any other candidate
distribution.  To determine which distribution has the best performance,
we utilize KS statistics to quantify the agreement between the empirical
distribution and the fitting distributions.
Figures~\ref{Fig:RI:Fit:Par}(a) and \ref{Fig:RI:Fit:Par}(d) show the KS
statistics of the five candidate distributions with respect to the mean
recurrence time $\tau_Q$ for the two stocks. For stock 000001 the
$q$-exponential distribution outperforms the other distributions and
possesses the smallest KS statistics for all $\tau_Q$. For stock 900956
the $q$-exponential distribution is best when $\tau_Q \le 60$, and the
stretched exponential distribution is best when $\tau_Q >
60$. Figure~\ref{Fig:RI:Fit:Par} shows plots of the characteristic
parameters of the five candidate distributions with respect to the mean
interval $\tau_Q$ for (b) stock 000001 and (e) stock 900956. Note that
all the fitting parameters of the five candidate distributions are
independent of the mean recurrence time $\tau_Q$ and exhibit a
horizontal line. This indicates that the distributions of the recurrence
intervals are not influenced by the threshold $Q$ when $\tau_Q$ is in
the $[20, 100]$ range. Figures~\ref{Fig:RI:Fit:Par}(c) and
\ref{Fig:RI:Fit:Par}(f) show a plot of the fitting parameters
$\lambda_x$, defined as $\tau_Q \lambda$, and $q$ as a function of
$\tau_Q$. Note that the fluctuations of $\lambda_x$ are not wide when
$\tau_Q > 40$ for stock 000001 and in the whole range of $\tau_Q$ for
stock 900956. These results indicate a scaling behavior in the
volatility recurrence intervals under different thresholds.

\begin{figure}[htb]
\centering
\includegraphics[width=5cm]{Fig_PDF_Fit_KS_000001.eps}
\includegraphics[width=5cm]{Fig_PDF_Fit_Par_000001.eps}
\includegraphics[width=5cm]{Fig_PDF_Fit_Par_qE_000001.eps}
\includegraphics[width=5cm]{Fig_PDF_Fit_KS_900956.eps}
\includegraphics[width=5cm]{Fig_PDF_Fit_Par_900956.eps}
\includegraphics[width=5cm]{Fig_PDF_Fit_Par_qE_900956.eps}
\caption{\label{Fig:RI:Fit:Par} (color online). Plots of fitting results
  for two stocks. (a-c) Stock 000001. (d-f) Stock 900956. (a, d) Plots
  of the KS statistics with respect to the mean interval $\tau_Q$. The
  KS statistics is used to describe the agreement between the empirical
  distribution and the fitting distributions. (b, e) Comparison of the
  characteristic fitting parameters for the five candidate
  distributions. (c, f) Plots of the fitting parameters $\lambda_x =
  \tau_Q\lambda$ and $q$ with respect to $\tau_Q$. }
\end{figure}

Our goal is to find a distribution that can approximate the distribution
of recurrence time. The best candidate distribution will be the one
that provides accuracy of fit and ease in estimating the hazard function
$W(t|\Delta t)$. Note that the power-law distribution with an
exponential cutoff provides the largest body of KS statistics for all
the stocks, and that the three-parameter Weibull distribution has fewer
KS statistics but requires that too many parameters be estimated. These
two distributions are excluded from the candidate list.  Because the
$q$-exponential distribution has the smallest body of KS statistics for
75.7\% of the fits, we use it to capture the distribution of the
recurrence intervals.  Even when the fits are of bodies of KS statistics
that are not the smallest, the $q$-exponential still provides a good
approximation of the recurrence time distribution. The $q$-exponential
distribution is also highly useful because it allows the derivation of
the analytical formula of the hazard function $W(t|\Delta t)$ from the
$q$-exponential formula.

\subsection{More on distribution parameter behaviors}

For each stock and each $\tau_Q$ value we fit the corresponding scaled
recurrence time $x$ with the $q$-exponential distribution and estimate
the distribution parameters $q$ and $\lambda_x$. To quantitatively check
the scaling behaviors we verify whether the parameters $q$ and
$\lambda_x$ are independent of $\tau_Q$ for the same stock and whether
they are the same across different stocks for the same $\tau_Q$. For
each stock we linearly regress the fitting parameters $q$ and
$\lambda_x$ with respect to $\tau_Q$. Although for $q$ vs $\tau_Q$ we
find that the maximum absolute slope is 0.006 and the mean slope over
all stocks $(5.68 \pm 7.77)\times 10^{-4}$, we do not observe comparable
small slopes for all stocks when fitting $\lambda_x$ vs $\tau_Q$. We
find 945 stocks with absolute slopes $<0.006$, which is the maximum
absolute slope of $q$ vs $\tau_Q$. The maximum value of the slope over
all stocks is 0.439 and the mean slope $0.0068 \pm 0.0213$ for
$\lambda_x$ vs $\tau_Q$. These results suggest that the parameter $q$ is
independent of the mean recurrence time $\tau_Q$ for all stocks and the
parameter $\lambda_x$ is independent of $\tau_q$ for half of the stocks,
which suggests that the scaling behavior in the recurrence intervals for
different $\tau_Q$ values exists in half of the stocks in the Chinese
markets ($945/1891 \approx 50\%$).

\begin{figure}[htb]
\centering
\includegraphics[width=5cm]{Fig_Fit_q_stock.eps}
\includegraphics[width=5cm]{Fig_Fit_q_Freq.eps}
\includegraphics[width=5cm]{Fig_Fit_lambda_Freq.eps}
\caption{\label{Fig:RI:Freq:FitPar} (color online). Plots of the
  estimated $q$ and $\lambda_x$ of $q$-exponential distribution for all
  stocks. (a) Plots of the fitting parameter $q$ for different
  stocks. There are five areas, three white areas are separated by two
  dark areas. For left to right, the first area represents the A-shares
  in Shenzhen market, the second area is the B-shares in Shenzhen
  market, the third area stands for ChiNext shares in Shenzhen market,
  the fourth area corresponds to the A-shares in Shanghai market, and
  the fifth area is the B-shares in Shanghai market. (b) Frequency of
  the fitting parameter $q$. (c) Frequency of the fitting parameter
  $\lambda_x$ for different values of $\tau_Q$.}
\end{figure}

Since the value of $q$ does not depend on $\tau_Q$, we average the
estimated $q$ for different $\tau_Q$ values for each stock and plot the
mean value of $q$ in Fig.~\ref{Fig:RI:Freq:FitPar}(a). Note the three
white areas separated by two shadow areas. The five areas from left to
right represent A-shares in the Shenzhen market, B-shares in the
Shenzhen market, ChiNext shares in the Shenzhen market, A-shares in the
Shanghai market, and B-shares in the Shanghai market. Note that there
are two groups of stocks, one with a relative smaller value of $\langle
q \rangle$ and the other with a much larger value of $\langle q
\rangle$. The group of stocks with the smaller $\langle q \rangle$ have
been traded in the market less than three years. The panel in
Fig.~\ref{Fig:RI:Freq:FitPar}(b) shows a plot of the frequency of
$\langle q \rangle$ in which there is a significant peak at $\langle q
\rangle\approx 1.3$, which corresponds to the average value of $q$ for
the stocks in the large $q$ group. Note also that there is a small peak
at $q\approx 1.1$, which is the mean value of $q$ for the stocks in the
small $q$ group. Because the fitting parameter $\lambda_x$ of half of
the stocks in our sample depends on $\tau_Q$, we plot the frequency of
$\lambda_x$ of all stocks for different values of $\tau_Q$ [see
  Fig.~\ref{Fig:RI:Freq:FitPar}(c)], instead of averaging the
$\lambda_x$ of different $\tau_Q$ for each stock. Note that the
distribution curve of $\lambda_x$ displays the same pattern across
different values of $\tau_Q$, the only difference being that the
distribution spanning range becomes wider when $\tau_Q$ increases. This
further indicates that the fitting parameter $\lambda_x$ is affected by
the mean recurrence time $\tau_Q$.

In order to determine whether the $q$-exponential parameters are
influenced by market state, i.e., bull or bear, we use a moving window
analysis to track the evolution of fitting parameters $q$ and
$\lambda_x$. Because stocks with trading records shorter than three
years have smaller $q$ values, we fix the window size at 48 months and
exclude stocks with trading periods shorter than 89 months.  We also
discard first-month trading data from the stocks remaining because
first-month records tend to be partial (i.e., do not span an entire
month) and the volatilities for new IPOs excessively large. We arrive at
1395 stocks as subject for our rolling window analysis.

\begin{figure}[htb]
\centering
\includegraphics[width=8cm]{Fig_Fit_q_MW.eps}
\includegraphics[width=8cm]{Fig_Fit_lambda_MW.eps}
\caption{\label{Fig:RI:Freq:FitPar:MW} (color online). Plots of the
  evolution of the estimated $q$ and $\lambda_x$ of $q$-exponential
  distribution for four chosen stocks. From the bottom panel to the up
  panel, the stock codes are 000001, 200002, 600220, and 900956,
  respectively. (a) Evolution of $q$. The dash line represents the mean
  value of $q$ across all windows. The solid line corresponds to the $q$
  of the entire period. (b) Evolution of $\lambda_x$ for different
  values of $\tau_Q$.}
\end{figure}

For each stock, we perform the same analysis in each moving window of
the series that we perform on the whole series. We first fit the
recurrence intervals for different values of $\tau_Q$ with the
$q$-exponential distribution and estimate the corresponding distribution
parameters $q$ and $\lambda_x$ in each window. We again find that the
parameter $q$ is independent of $\tau_Q$ in each window for each stock,
but the parameter $\lambda_x$ does not exhibit this behavior. We average
the $q$ for different values of $\tau_Q$ in each window and plot the
mean $q$ as a function of the time, which corresponds to the last month
of each moving window [see Fig.~\ref{Fig:RI:Freq:FitPar:MW}(a)]. For the
sake of comparison, the estimated $q$ of the entire period is a solid
line and the mean value of $q$ across all the windows is a dashed
line. Note that for the four stocks there is a big gap between the solid
and dashed line. Note also that the $\langle q \rangle$ curves also
exhibit significant fluctuations along the time
axis. Figure~\ref{Fig:RI:Freq:FitPar:MW}(b) shows the evolution of
$\lambda_x$ for different values of $\tau_Q$. Note that the trajectory
of $\lambda_x$ exhibits the same trend as the trajectory of $\langle q
\rangle$. It is not clear whether these trends reflect stock price
tends, but our results demonstrate that the estimated parameters of the
$q$-distribution are influenced by market status and thus may be treated
as a individual risk factor when explaining market returns.

\section{Results of the large volatility prediction}
\label{S1:Risk}

\subsection{Hazard probability}

The recurrence intervals between volatilities exceeding a threshold $Q$
are well approximated by the $q$-exponential distribution, and that
gives us the formula of interval distribution $p(t) = (2-q) \lambda
[1+(q-1) \lambda t]^{-\frac{1}{q-1}}$. By substituting this equation
into equation~(\ref{Eq:Wq}), we obtain
\begin{equation}
 W_q(\Delta{t}|t)= 1- \left[ 1 + \frac{(q-1)\lambda \Delta
     t}{1+(q-1)\lambda t} \right]^{1-\frac{1}{q-1}}.
 \label{Eq:Wq:qe}
\end{equation}

If we designate the top 1\% of volatility values (corresponding to the
mean recurrence time $\tau_Q = 100$) to be extreme events, we can
estimate the hazard probability $W_q(\Delta{t},t)$ in $t$ when fixing
$\Delta t$. Figure~\ref{Fig:RI:Wdtt} shows the estimated hazard
probability for two stocks when $\Delta t = 1$, 5, and 10. Note that
both the solid lines indicating the analytical solution
equation~(\ref{Eq:Wq:qe}) and the markers indicating the empirical data
decrease slowly and in each panel are in good agreement. The decreasing
trend of $W(\Delta t |t)$ is consistent with the clustering behavior in
the volatility series. Note that the hazard probability $W(\Delta t |t)$
is universal and can be used to estimate the risk in any kind of time
series.

\begin{figure}[htb]
\centering
\includegraphics[width=5cm]{Fig_Wdtt_000001_1.eps}
\includegraphics[width=5cm]{Fig_Wdtt_000001_5.eps}
\includegraphics[width=5cm]{Fig_Wdtt_000001_10.eps}
\includegraphics[width=5cm]{Fig_Wdtt_900956_1.eps}
\includegraphics[width=5cm]{Fig_Wdtt_900956_5.eps}
\includegraphics[width=5cm]{Fig_Wdtt_900956_10.eps}
\caption{\label{Fig:RI:Wdtt} (color online). Plot of hazard probability
  $W(\Delta t|t)$ for two stocks and $\Delta t = 1$, $5$, and
  $10$. (a-c) Stock 000001. (d-f) Stock 900956.}
\end{figure}

\subsection{Predicting large volatilities}

To forecase large volatility events in a volatililty series we first
calculate the recurrence time for a given $\tau_Q$, e.g., $\tau_Q =
100$, and estimate the distribution of parameters. The hazard
probability $W(1|t)$ that an extreme event will occur in the next period
is determined using equation~(\ref{Eq:Wq:qe}) and the estimated distribution
parameters. Once the hazard probability breaks though the predefined
threshold $Q_p$, an alarm will be triggered, warning that a large
volatility event is immanent. Figure~\ref{Fig:RI:PreExtE}(a) plots a
subseries of the volatility values and highlights events above threshold
$Q$ in the top panel, which correspond to the mean recurrence time
$\tau_Q$. The risk probability $W(1|t)$ is shown in the bottom
panel. Note that $W(1|t)$ decreases as time $t$ elapsed from the last
large volatility event increases.  Threshold $Q_p$ is plotted as a
horizontal line to show the activating alarm process. By varying $Q_p$
within a $[0, 1]$ range, we obtain all pairs of $(A, D)$.

Figure~\ref{Fig:RI:PreExtE}(b) shows the ROC curves for ten stocks. Note
that all ten curves are above the dashed line $D=A$, indicating that our
prediction is not random.  The ten curves do not overlap, indicating
that the accuracy of this prediction algorithm varies for different
stocks. For the same false alarm $A$ (where a vertical dashed line is
plotted at $A=0.1$ as an example), stock 900901 has the highest correct
prediction rate and stock 300066 the lowest. We also calculate the
correct prediction rate $D_{A=0.1}$ for all stocks at the false alarm
level of 0.1. Figure~\ref{Fig:RI:PreExtE}(c) shows the frequency plots
of $D_{A={\rm{0.1}}}$. Note that the peak is centered at $\approx
0.2$, 10\% higher than a random prediction. We also find there are
nine stocks with $D_{A=0.1} > 0.4$. The top three values are 0.7, 0.68,
and 0.60 for stocks 000529, 000557, and 000592, indicating that our
forecasting algorithm can accurately predict the large volatilities of
these three stocks. Previous research has indicated that the efficiency
of the algorithm is primarily influenced by the linear and nonlinear
memory behavior in the original volatilities
\citep{Bogachev-Bunde-2011-PA}. Our results here indicate that the
behavior of stocks with stronger memory behaviors, such as volatility
clustering and multifactality, could be more accurately predicted using
our algorithm. Our algorithm only takes into consideration the
probability distribution of recurrence intervals, but if the memory
behavior of recurrence intervals were also included we believe that its
predictive accuracy would be greatly inhanced.

\begin{figure}[htb]
\centering
\includegraphics[width=5.5cm]{Fig_Volatility_HazardProbability_600220.eps}
\includegraphics[width=5cm]{Fig_RI_PreEve_AD.eps}
\includegraphics[width=5cm]{Fig_RI_PreEve_DisD.eps}
\caption{\label{Fig:RI:PreExtE} (color online). Prediction of large
  volatilities. (a) Plots of a representative volatility series in the
  top panels and hazard probability $W(1|\Delta t)$ in the bottom
  panels. (b) Plots of ROC curves for ten stocks. (c) Distribution plots
  of the correct prediction rates at the false alarm level of 0.1 for
  all stocks.}
\end{figure}

\section{Conclusion}
\label{S1:Conclusion}

In this work, we have utilized a decision-making algorithm to forecast
the occurrence of large volatilities in Chinese stock markets based on
the hazard probability, which is derived from the distribution of
recurrence intervals between the volatilities exceeding a threshold
$Q$. By fitting the volatility recurrence intervals by means of five
candidate distributions and comparing their KS statistics, we have found
that the volatility recurrence intervals are well approximated by the
$q$-exponential distribution. The fitting parameter $q$ is found to be
independent of the mean recurrence time $\tau_q$, which is at a
one-to-one correspondence with the threshold $Q$, for all the stocks in
our sample. However the parameter $\lambda_x$ of the $q$-distribution
does not exhibit the same behavior as $q$. For half of the stocks,
$\lambda_x$ exhibits a strong dependence on $\tau_Q$. Using a moving
window analysis, we have found that both parameters are influenced by
market status and exhibit the same trend with the evolution of
time. This behavior may have potential applications for explaining stock
return volitility.

Using the $q$-distribution formula, we have derived an analytical
solution of the hazard probability $W(\Delta t|t)$ of the next large
volatility event above the threshold $Q$ within a short time interval
$\Delta t$ after an elapsed time $t$ from the last large volatility
above $Q$. This analytical solution $W(\Delta t|t)$ is in good agreement
with the empirical risk probability derived from real stock data. We
have adopted a decision-marking algorithm and have used the hazard
probability to forecast large volatilities. At the false alarm level of
0.1, we have found that the average correct predicting rate is 0.2 for
all stocks. We have also found that there are three stocks with a
correct predicting rate is greater than 0.6, indicating that our
predicting algorithm is accurate in forecasting the large
volatilities. Our findings may shed new light on our understanding of
extreme volatility behavior and may have potential applications in
managing stock market risk.

\section*{Acknowledgements}

Z.-Q.J. and W.-X.Z. acknowledge support from the National Natural
Science Foundation of China (71131007), Shanghai ``Chen Guang'' Project
(2012CG34), Program for Changjiang Scholars and Innovative Research Team
in University (IRT1028), China Scholarship Council (201406745014) and
the Fundamental Research Funds for the Central
Universities. A.C. acknowledges the support from Brazilian agencies
FAPEAL (PPP 20110902-011-0025-0069/60030-733/2011) and CNPq (PDE
20736012014-6).

\bibliographystyle{elsarticle-num}


\section{Introduction}
\label{sec:intro}

Fingering convection (also called thermohaline convection) was first discussed in the astrophysical context by \citet{Ulrich1972} \citep[see also][]{Ulrich1971}. The basic fingering instability, which occurs in systems that are compositionally unstably stratified and thermally stably stratified, is common in stars whose surface is polluted by high mean molecular weight material. This can happen through the accretion of material from infalling planets or debris disks \citep{vauclair2004mfa,garaud11,Dealal2013} or from more evolved companion stars \citep{stancliffe2007cem,Angelou2012}. Fingering convection can also result from the radiative levitation of heavy elements such as iron from deeper regions upwards in relatively high-mass stars \citep{theado09,Zemskovaal14}, or from off-center nuclear burning in post main sequence stars \citep{CharbonnelZahn07,denissenkov2010,DenissenkovMerryfield2011}.

In the typical conditions encountered in stellar interiors, instability to fingering convection only depends on the value of the density ratio $R_0$, defined as
\begin{equation}
R_0 = \frac{\alpha \left( \frac{dT_0}{dr} - \frac{dT_{\rm ad}}{dr} \right) }{\beta \frac{d\mu_0}{dr} } \mbox{   ,   }
\end{equation}
where $T_0(r)$ is the local temperature profile in the stellar region considered, $T_{\rm ad}(r)$ is the temperature profile that region would have were it to be adiabatically stratified, $\mu_0(r)$ is the local mean molecular weight profile, and $\alpha$ and $\beta$ are derivatives of the equation of state defined as
\begin{equation}
\alpha = - \frac{1}{\rho} \left(\frac{\partial \rho}{\partial T} \right)_{P,\mu} \mbox{   ,   } \beta = \frac{1}{\rho}\left( \frac{\partial \rho}{\partial \mu}\right)_{P,T}  \mbox{   .  }
\end{equation}
The density ratio is therefore the ratio of the density gradient due to the stabilizing thermal stratification to the density gradient due to the destabilizing compositional stratification.
 Note that $R_0$ is also commonly written as
\begin{equation}
R_0 = \frac{\delta  (\nabla - \nabla_{\rm ad})}{\phi \nabla_\mu} \,
\end{equation}
where $\nabla$, $\nabla_{\rm ad}$ and $\nabla_{\mu}$ have their usual astrophysical definitions, and
\begin{equation}
\delta = - \left( \frac{\partial \ln \rho}{\partial \ln T} \right)_{P,\mu} =  \alpha T \mbox{  and  } \phi =  \left( \frac{\partial \ln \rho}{\partial \ln \mu} \right)_{P,T}  = \beta \mu  \mbox{   .  }
\end{equation}

As shown by \citet{bainesgill1969}, linear instability to fingering convection only occurs when
\begin{equation}
1 < R_0 < R_{\rm crit} = \tau^{-1} = \frac{\kappa_T}{\kappa_\mu} ,
\end{equation}
where $\kappa_T$ and $\kappa_\mu$ are the thermal and compositional diffusivities respectively.  The lower limit ($R_0 = 1$) corresponds to the onset of standard overturning convection, while the upper limit ($R_0 = R_{\rm crit}$) corresponds to the critical point of marginal stability to fingering convection.  As shown by \citet{Medranoal15}, $\tau$ varies between $10^{-9}$ and $10^{-5}$ in the interiors of main sequence stars, and increases only up to about $10^{-3}$ in degenerate regions of more evolved stars (e.g. red giants or white dwarfs). This shows that $R_{\rm crit}$ is always much greater than one, so that stellar interiors can very easily be destabilized: even a tiny adverse compositional gradient can trigger fingering instabilities. As a result, fingering convection is very common, and should be taken into account when considering mixing in stellar evolution models.

An important difference in the behavior of fingering instabilities in astrophysical and geophysical systems comes from the value of the Prandtl number (${\rm Pr} = \nu/\kappa_T$, where $\nu$ is the kinematic viscosity): this parameter is usually larger than one in most geophysical fluids, but is asymptotically small \citep[typically, at most one order of magnitude larger than $\tau$; see][]{Medranoal15} in stellar interiors. As a result, fingering structures which are typically fairly laminar in the geophysical context are instead very turbulent in the astrophysical one. This complicates the study of the saturation of the instability, which is necessary to estimate the rate of heat and compositional transport induced by fingering convection.

In the past few years, thanks to the development of fast parallel algorithms and advances in supercomputing, it has become possible to follow the nonlinear development of fingering instabilities in three dimensional (3D) simulations at relatively low diffusivity ratio and low Prandtl numbers, albeit not yet at actual stellar values of these parameters. The first of such studies were by \citet{denissenkov2010} using 2D simulations and by \citet{Traxleral2011} and \citet{Brownal2013} in 3D. The 3D simulations covered a wide range of parameter space with ${\rm Pr}$ and $\tau$ varying from 0.3 down to 0.01, and the density ratio was varied across the instability range (from 1 to $\tau^{-1}$). Using them, \citet{Brownal2013} were able to propose an analytical model for turbulent transport by fingering convection at low ${\rm Pr}$ and $\tau$. This model has only one free parameter, which can be calibrated using the 3D simulations.
Once calibrated, the \citet{Brownal2013} model correctly predicts the transport rate for both heat and composition for {\it all} available simulations in which $\tau < {\rm Pr}$ within a factor of 2 at worst, often much better.

Whether this model remains as accurate for $\Pr$ and $\tau$ lower than $10^{-2}$ still remains to be determined. However, cubic-domain 3D simulations at very low ${\rm Pr}$ and $\tau$ rapidly become computationally prohibitive since the resolution must  be increased as ${\rm Pr}$ and $\tau$ decrease, in order to fully resolve the various boundary layers that form at the edge of the fingers. In addition, and as discussed by \citet{Traxleral2011b} and \citet{Medranoal15}, secondary large-scale instabilities can develop spontaneously from fingering convection, and the latter substantially affect the turbulent transport properties of the system. To study them requires very large computational domains, containing at the very least 20 to 30 wavelengths of the fastest-growing fingering mode in both horizontal and vertical directions, and often much more than that. Again, using such large domains in 3D, together with low ${\rm Pr}$ and $\tau$, is currently computationally prohibitive.

For both reasons, it is very tempting to go back to 2D simulations either for very low diffusivity studies, or for very large domain studies. In fact, 2D fingering convection simulations used to be common in the oceanographic literature only 20 years ago \citep[e.g.][]{shen1995,stern2001sfu}, and are still often in use today \citep{SimeonovStern2007,radko2008ddm,Sreenivasal09,SinghSrinivasan14} in that context. \citet{stern2001sfu} showed by comparing 2D and 3D simulations at ${\rm Pr} = 7$ that the typical turbulent fluxes estimated from 2D simulations do not vary from those obtained in 3D simulations by more than a factor of a few. \citet{DenissenkovMerryfield2011} claim to have run both 2D and 3D simulations of fingering convection for ${\rm Pr}$ and $\tau$ both $O(10^{-6})$, and to have found that the effective turbulent diffusivities of composition are very similar in the two cases.

However, the work of \citet{Radko2010} casts strong doubts on the claims of \citet{DenissenkovMerryfield2011}. Indeed, he showed both analytically and numerically that the saturation of the fingering instability at low ${\rm Pr}$ in 2D can be attributed to the development of strong shear layers that destroy the fingers and completely change the overall behavior of the induced turbulence. A similar process was not seen in the 3D simulations of  \citet{RadkoSmith2012} and \citet{Brownal2013}. This raises a number of practical questions: (1) Is shear indeed present in low Prandtl number fingering \citep[as claimed by][]{Radko2010}, or not \citep[as implied by][]{DenissenkovMerryfield2011}? (2) More generally, under which circumstances, if any, can 2D simulations be used to study fingering convection at low Prandtl number? (3) If 2D simulations are not appropriate, then how thick does a domain have to be in order to approximate 3D results appropriately? In this work, we answer all three questions.

Section \ref{sec:model} outlines our model equations, boundary conditions, and domain geometry. In Section \ref{sec:highR0}, we first compare 2D and 3D results in high density ratio simulations. In Section \ref{sec:lowR0}, we present a similar comparison but in low density ratio simulations. As we shall demonstrate, in both cases 2D and 3D results are fundamentally different, but for different reasons. Section \ref{sec:ccl} discusses our findings and proposes acceptable compromises in terms of running 3D fingering simulations as cheaply as possible.

\section{Model}
\label{sec:model}

In what follows, we study the dynamics of homogeneous fingering convection in the Boussinesq approximation, as introduced for the purpose of numerical simulations by \citet{shen1995} \citep[see also][]{Radko2003mlf}. In this formalism, the background temperature and composition profiles are assumed to be linear, with locally constant gradients $dT_0/dr$ and $d\mu_0/dr$ (see Section \ref{sec:intro}), while all perturbations to that background (including the temperature $T$, the pressure $p$, the velocity field $\vec{u}=(u,v,w)$ and the mean molecular weight $\mu$) are assumed to be triply periodic in the domain considered.

The non-dimensional governing equations are given by:
\begin{eqnarray}
&& \frac{1}{\rm{Pr}}\left(\frac{\partial \vec{u}}{\partial t} + \vec{u}\cdot\nabla \vec{u}\right) = -\nabla p + (T-\mu) \vec{e}_z + \nabla^2 \vec{u}, \nonumber \\
&& \frac{\partial T}{\partial t} + \vec{u}\cdot\nabla T + w  = \nabla^2 T,\nonumber \\
&& \frac{\partial \mu}{\partial t} + \vec{u}\cdot\nabla \mu + \frac{1}{R_0} w  = \tau \nabla^2 \mu,\nonumber \\
&& \nabla \cdot \vec{u} = 0
\label{eq:nondim}
\end{eqnarray}
where ${\rm Pr}$, $\tau$, and $R_0$ were introduced in Section \ref{sec:intro}.
 To arrive at these equations, we have used the following units:
\begin{eqnarray}
&&  \left[l\right] = d = \left(\frac{\kappa_T \nu} {\alpha g\left|\frac{dT_0}{dr} - \frac{dT_{\rm ad}}{dr}\right|  }\right)^{1/4} \mbox{ as the unit distance} , \nonumber \\
&& \left[t\right] = \frac{d^2}{\kappa_T} , \quad \, \left[u\right] = \frac{\kappa_T}{d} \mbox{ as the unit time and velocity respectively} , \nonumber \\
&& \left[T\right] =  \left|\frac{dT_0}{dr} - \frac{dT_{\rm ad}}{dr}\right| d \mbox{ as the unit temperature} , \quad \, \left[\mu\right] = \frac{\alpha }{\beta}  \left[T\right]  \mbox{ as the unit composition}
\end{eqnarray}
where $g$ is the local gravity, and all other quantities were defined in Section \ref{sec:intro}.

In what follows we consider two possible types of domains: 2D domains in the $(x,z)$ plane of size $(250d \times 250d)$, and 3D domains of size $(250d \times L_y \times 250d)$, where $L_y$ shall be varied. We use a fairly large size in the $(x,z)$ plane to allow for the development of large-scale structures in the simulations, should they naturally occur.
We set the parameters ${\rm Pr}$ and $\tau$ to be both equal to $0.03$, well below unity. The density ratio, on the other hand, is varied to detect and study the existence of different dynamical regimes across the fingering instability range as appropriate. The code used to solve the set of equations (\ref{eq:nondim}) under triply-periodic boundary condition is the pseudo-spectral code described in \citet{Traxleral2011b}. In all cases, we use a resolution of $256$ spectral modes per $250d$ (corresponding to an effective resolution of about 3 grid points per $d$). The same resolution is used in all directions.

\section{Comparison of 2D and 3D runs at high density ratio}
\label{sec:highR0}

We begin our investigation by considering a ``high'' density ratio case, with $R_0 = 5$. While still significantly below $R_{\rm crit} = 33.3$ (at these parameters), this value is already large enough to be in the sheared regime described by \citet{Radko2010} (see Section \ref{sec:intro}), but small enough to ensure that resolving the buoyancy frequency (which is equal to $\sqrt{\Pr (R_0 - 1)}$ in these units) remains computationally manageable. Indeed, gravity waves are present in this system, and the time step must be small enough to fully resolve their oscillation period.

\subsection{Results of 2D simulations}
\label{sec:2Dhi}

Figure \ref{fig:2DhiR0snaps} shows the evolution of the mean molecular weight perturbation $\mu$ and of the horizontal component of the velocity field $u$ in the 2D run, from very early times to much later times. We see, successively, the development of the basic fingering instability (top), its saturation (middle), and finally, the emergence of strong alternating shear layers and their effect on the fingers (bottom). As seen in Figure \ref{fig:layermerger}, the shear layers evolve slowly with time: weak layers appear to be pushed closer to one another, and eventually merge into stronger ones. As discussed in Section 1, the spontaneous emergence of large-scale shear layers in 2D fingering convection at low Prandtl number had already been discovered and studied by \citet{Radko2010}. We therefore confirm his findings here.

\begin{figure}
\centerline{\includegraphics[width=0.8\textwidth]{Figure1.pdf}}
\caption{Snapshots of the compositional perturbations (left) and of the horizontal velocity in the $x$ direction (right) at three times (see the evolution also shown in Figure \ref{fig:Nuurms}): near the onset of the fingering instability ($t = 280$, top), early after its saturation ($t = 530$, middle) and at later time once the shear layers have begun to develop ($t =2500$ bottom). Note how the horizontal shear distorts the fingers. }
\label{fig:2DhiR0snaps}
\end{figure}
\begin{figure}
\centerline{\includegraphics[width=0.8\textwidth]{Figure2.pdf}}
\caption{Horizontally-averaged horizontal velocity profile in the 2D simulation of fingering at $R_0 = 5$ as a function of time. The layers rapidly merge down to four, after which a quasi-steady periodic state appears to be reached. The oscillations in the intensity of the shear are clearly visible.}
\label{fig:layermerger}
\end{figure}

Interestingly, we further find that the nonlinear interaction of the fingers and the shear drive what appears to be relaxation-oscillations. The latter are most noticeable in Figure \ref{fig:Nuurms}, when viewed in terms of the total compositional flux $\langle w\mu \rangle$ and of the r.m.s. horizontal velocity $u_{\rm rms} = \langle u^2 \rangle^{1/2}$ (where the brackets $\langle \cdot \rangle$ denotes a volume average over the entire domain). The oscillations appear around $t = 5000$, then grow in amplitude as the various shear layers merge. By $t = 15000$, the oscillation pattern becomes very clear: efficient fingering (characterized by a large vertical flux) drives horizontal shear (noticeable in the increase of $u_{\rm rms}$), which first distorts the fingers then eventually completely suppresses fingering convection (hence the drop in the turbulent flux). When this happens, the mechanism that drives the shear dies out, and the shear gradually disappears. The fingers are eventually allowed to grow again, and the cycle repeats. Note that the amplitude and period of the cycle appear to depend on the number of shear layers present. As the simulation progresses and the number of layers decreases through mergers, the cycle lengthens and its amplitude grows, until a quasi-periodic state is reached.
\begin{figure}
\centerline{\includegraphics[width=0.5\textwidth]{Figure3.pdf}}
\caption{Compositional flux and r.m.s. horizontal velocity in the 2D simulation of fingering at $R_0 = 5$, illustrating the mechanism responsible for the relaxation oscillations.}
\label{fig:Nuurms}
\end{figure}

\subsection{Results of 3D simulations of various $L_y$}
\label{sec:3Dhi}

Even though these fascinating dynamics appear in 2D simulations of fingering convection, it is easy to show that the latter are spurious and do not exist in computational domains that are ``sufficiently 3D". Figure \ref{fig:2D3DhiR0} shows the turbulent compositional flux $\langle w\mu \rangle$ and the r.m.s. horizontal velocity $u_{\rm rms}$ as a function of time for simulations that are perfomed at exactly the same parameters and have $R_0 = 5$, as in Section \ref{sec:2Dhi}, but now in 3D with different domain thicknesses $L_y$. The early stages of the 2D simulation discussed in the previous section are included for comparison.

It is quite clear that the simulation for which $L_y = 4d$ still behaves as if it were 2D, that is, with a turbulent flux that rapidly decreases with time after the initial saturation, and a corresponding r.m.s. horizontal velocity that increases with the generation of shear layers.  By contrast, simulations with $L_y = 15d$ and larger all behave in quantitatively similar ways, achieving the same well-defined quasi-steady state post-saturation, with substantially higher flux and lower horizontal velocities than in 2D and no obvious quasi-periodic oscillations. The run with $L_y = 8d$ appears to have the same character as the wider 3D domains but somewhat underestimates the flux.

The significance of the  cutoff-sizes $L_y = 8d$ and $L_y = 15d$ becomes clearer if one notes that, at the parameter values selected, the typical wavelength of the fastest-growing mode of fingering convection is about 8.5d (one wavelength containing one up-going and one down-going finger). We conclude from this experiment that at low Prandtl number and high density ratio: (1) a domain has to contain at least 1 wavelength of the fastest-growing fingering mode in the third dimension to exhibit dynamics that are fully 3D instead of the spurious shear layers and associated relaxation oscillations that are observed to exist in 2D; (2) a domain has to contain at least 2 wavelengths of the fastest-growing mode to yield quantitatively accurate results for the turbulent fluxes and turbulent intensity of fingering convection. These results are generic:  similar quantitative findings have been obtained for other values of ${\rm Pr}$, $\tau \ll 1$ and for sufficiently large $R_0$.

\begin{figure}
\centerline{\includegraphics[width=\textwidth]{Figure4.pdf}}
\caption{Compositional flux (left) and r.m.s. horizontal velocity (right) for fingering simulations at $R_0 = 5$, of varying domain thicknesses (see legend for detail). Runs with very thin domain behave as if they were 2D, while thick enough domains behave as in a 3D manner that rapidly becomes independent of the domain size. }
\label{fig:2D3DhiR0}
\end{figure}

\subsection{Can the artificial shear be prevented? }

The results of Figure \ref{fig:2D3DhiR0} strongly suggest that the substantial horizontal shear and the associated quenching of the vertical fluxes observed in the 2D case are artificial. However, given the vast savings in computational time between 2D and 3D simulations, it is tempting to wonder if one could still recover acceptable 2D solutions (i.e. solutions whose dynamical behavior and turbulent transport rate are comparable to the 3D solutions) by suppressing the shear. Figures \ref{fig:nic1} and \ref{fig:nic2} demonstrate that this does not work.  Since the code we are using is spectral, it is easy to zero out the mean component of the horizontal flow at every time step, effectively suppressing the shear whilst allowing all other modes to evolve normally.  We performed this experiment on the 2D case at $R_0=5$ discussed in Section \ref{sec:2Dhi}. We compare below these new results to those originally obtained without shear-suppression, and to the narrow 3D case with $L_y = 15d$ discussed in Section \ref{sec:3Dhi} that appeared representative of the true solution.  Finally, in order to ensure that shear-suppression does not artificially affect the dynamics of fingering convection in 3D, we ran a final experiment suppressing the shear in the same manner on another wider 3D case (with domain size $250d\times50d\times250d$) at the same parameters.

Figure \ref{fig:nic1} shows the thermal flux $\langle wT \rangle$ for each of these cases.   As seen before, the narrow 3D case settles down to a much higher flux than the 2D case, which exhibits relaxation oscillations.  The wider 3D case with the shear suppressed has the same flux as the narrower 3D one, showing that shear does not play a significant role in the 3D simulations.  The 2D case with the shear suppressed, however, does something altogether different from all the other cases.  Although it initially settles into a mode with roughly the right mean flux, the fluctuations about the mean are much larger than in the 3D cases.  Furthermore, at around $t=40000$, the 2D shear-suppressed case jumps to a solution with a much higher flux.

This rather peculiar behavior is best understood by looking at snapshots of the simulations. Figure \ref{fig:nic2}, which shows slices of the horizontal velocity field, demonstrates the physical differences between these various cases.  In 3D, the solution consists of a homogeneous field of small-scale fingers, whether the shear is removed or not.  In the original 2D case, the presence of strong shear discussed in Section \ref{sec:2Dhi} is obvious.  In 2D with the shear suppressed, it appears that the solution attempts to create a velocity pattern that is as close to the regular 2D case (with strong mean shear) as possible, given the model constraints. After a short transient phase, strong bands of flow appear. They are not exactly horizontal, as that is disallowed, but are slightly angled upwards instead.  To begin with, the specific angles are somewhat random over space. The corresponding vertical flux is significantly larger than in the regular 2D case, and somewhat comparable to the one observed in 3D.  Ultimately though, the flow self-organizes at a single well-defined angle, and acts as an extremely efficient periodic ``escalator" that significantly enhances the vertical transport. Despite being quite interesting, however, these flows remain clearly artificial compared to the 3D cases. We are therefore forced to conclude that there is no easy fix to the problem, and that 3D simulations are the only way to get reliable results.

\begin{figure}[h]
\centerline{\includegraphics[width=0.8\textwidth]{Figure5.pdf}}
\caption{Comparison of time evolution of thermal fluxes for 2D and 3D simulations at $R_0=5$ with and without mean shear flows suppressed.  Fluxes are shown for a purely 2D simulation, a purely 2D simulation with shear suppressed, a narrow 3D simulation that contains 2 finger wavelengths in the narrow direction, and a wider 3D simulation with the shear suppressed.  Clearly, whilst both 3D simulations agree well, both 2D simulations do not reflect the same behavior at all.}
\label{fig:nic1}
\end{figure}
\begin{figure}[h]
\centerline{\includegraphics[width=\textwidth]{Figure6.pdf}}
\caption{Slices in the $(x,z)$ plane of the horizontal velocity for the simulations shown in Fig. \ref{fig:nic1} at late times in each simulation.  In 2D, strong shear flows are seen, and in 2D with mean shear suppressed, an ``elevator'' solution is observed that is very close to the sheared 2D case but which is periodic.  None of this behavior is seen in 3D, where the suppression of shear makes no difference since significant shear is not present.}
\label{fig:nic2}
\end{figure}

\section{Comparison of 2D and 3D runs at low density ratio}
\label{sec:lowR0}

We now turn to low density ratio systems, and illustrate our findings with a case which has $R_0 = 1/0.9 \simeq 1.11$, and ${\rm Pr} = \tau = 0.03$ as above.
Figure \ref{fig:2D3DloR0snaps} shows snapshots of a 2D simulation, and of a 3D simulation done in the domain of thickness $L_y = 125d$, both taken around $t = 500$. Again, we see a notable difference between the 2D and the 3D case, but the nature of this difference is clearly not the same as the one discussed in Section \ref{sec:highR0}.
In the 2D case, we see large plumes that are reminiscent of oversized turbulent fingers. It is worth remembering that the typical basic fingering mode width is about 1/30 of the domain size, so the structures seen here (which have a size of about 1/3 of the domain) are much larger than basic fingers. The 3D case, by contrast, exhibits very different dynamics. We note a clear separation of scales between small-scale fingers (which have more or less the same size as the basic instability) and domain-scale gravity waves that modulate the fingering field.

The emergence of these large-scale gravity waves is well-understood, and can be attributed to the collective instability first discussed by \citet{stern1969cis} in the oceanographic context, later found in simulations of low-Prandtl number fingering convection by \citet{Brownal2013}, and systematically studied by \citet{Medranoal15}. The collective instability is a mean-field instability, excited through a positive feedback loop between the large-scale gradients of temperature and composition, and the respective turbulent fluxes induced by these large-scale gradients\footnote{As discussed by \citet{Radko13}, the collective instability can be interpreted as the turbulent analog of Oscillatory Double-Diffusive Convection, a linear instability that normally takes place in semiconvective regions \citep{Kato66}.}. These large-scale gravity waves do not appear to be excited in 2D, which is somewhat surprising since the collective instability theory does not a priori need to be 3D to operate.

\begin{figure}
\centerline{\includegraphics[width=\textwidth]{Figure7.pdf}}
\caption{Snapshots of the compositional perturbation for a simulation with $R_0 = 1.11$, in 2D (left) and in a thick 3D domain (right). Note the presence of large-scale gravity waves in the 3D case, absent in the 2D case.  Also note the significant difference in the maximum amplitude of the compositional perturbations in each case.}
\label{fig:2D3DloR0snaps}
\end{figure}

A more quantitative way of seeing the scale separation inherent to the collective instability in 3D, and its absence in 2D, is by inspection of the kinetic energy spectrum. This is shown in Figure \ref{fig:specs}  as a function of the vertical wavenumber, for the same simulations and at the same times as the ones shown in Figure \ref{fig:2D3DloR0snaps}. The spectrum of the 2D run clearly peaks around $k  = 0.08 d^{-1}$, which corresponds to a typical wavelength of about $80d$, or in other words, about 10 times the wavelength of the fastest-growing fingering mode. This is more-or-less the scale of the plumes observed in the snapshot. Between that value of $k$ and the energy injection scale ($k \sim 1$) we observe that the kinetic energy spectrum is close to a power law with index $-2$, which may indicate the presence of an inverse energy cascade. The 3D run, by contrast, has a kinetic energy spectrum that peaks at the lowest possible wavenumber $k = 0.025d^{-1}$, which corresponds to a wavelength commensurate with the domain size. This is indeed the scale of the collective instability mode observed in the snapshot. The energy spectrum drops sharply between $k = 0.025d^{-1}$ and $k = 0.1d^{-1}$, then more gently between $k \sim 0.1d^{-1}$ and $k \sim d$ (the energy injection scale).  Beyond $k=d$, both 2D and 3D simulations have the same energy spectrum which is dominated by small-scale diffusive processes, showing that they are basically the same in 2D and 3D.

\begin{figure}
\centerline{\includegraphics[width=0.7\textwidth]{Figure8.pdf}}
\caption{Kinetic energy spectrum as a function of the vertical wavenumber $k$, for the two simulations shown in Figure \ref{fig:2D3DloR0snaps}. The linear instability injection scale is estimated from the horizontal wavenumber $l$ of the fastest-growing mode, and assuming that $ k \sim l$. }
\label{fig:specs}
\end{figure}

Finally, Figure \ref{fig:2D3DloR0} shows the turbulent compositional flux for simulations where the thickness of the domain $L_y$ is varied from $0$ to $125d$ as we did in Figure \ref{fig:2D3DhiR0} for the high density ratio case. This time, we see a clear difference in behavior between all the 3D simulations and the 2D simulation. The 2D case shows a continued growth of the efficiency of turbulent transport that only saturates quite late (after $t = 500$) at a fairly high mean value. All the 3D simulations, by contrast, show an early saturation of the primary fingering instability around $t = 150$, then a secondary growth phase associated with the growth of the collective instability, which later also saturates but at a much lower mean value than the 2D case.
By contrast with the high density ratio case described in the previous Section, even the $L_y = 4d$ simulation seems to be qualitatively more compatible with the full 3D case than the 2D case. However, in order to get quantitatively accurate flux measurements just after the saturation of the fingering instability but prior to the onset of the collective instability (i.e. roughly between $t = 150$ and $t = 250$), it is necessary to choose $L_y \ge 15d$, just as in the high density ratio case. In other words, simulations in domains that are at least 2 wavelengths of the fastest-growing mode (which is also about $8.5d$ at this value of $R_0$) are necessary to get adequate estimates of the fingering fluxes at low Prandtl number.

Of course, even if the fingering dynamics are correctly accounted for, a very narrow domain may no longer be adequate once larger-scale structures appear. In the $L_y = 15d$ domain, the large-scale gravity waves excited by the collective instability are artificially forced to be 2D, for instance. How this affects their saturation amplitude, and their overall behavior, remains to be determined. Similarly, in the work of \citet{Zemskovaal14}, the thin domain used severely restricts the dynamics of the convective layer that eventually forms, to the extent that the results {\it in the convective phase} are not necessarily reliable\footnote{see, for instance, the strong vortices that are clearly visible in their Figure 10, a classic sign of quasi-2D dynamics that doesn't persist in 3D}. In other words, while a thin domain is a good compromise to study homogeneous fingering convection, a full 3D domain with aspect ratio of order unity may remain necessary to properly model any large-scale dynamics that naturally arise.

\begin{figure}
\centerline{\includegraphics[width=0.8\textwidth]{Figure9.pdf}}
\caption{Compositional flux for $R_0 = 1.11$, for various domain thicknesses (see legend for detail). The 2D and 3D runs behave very differently. }
\label{fig:2D3DloR0}
\end{figure}

\section{Discussion and prospect}
\label{sec:ccl}

In this work we have compared 2D and 3D simulations of fingering convection at low Prandtl number to determine under which circumstance 2D models can be used to approximate  3D systems. We have found that 2D simulations always suffer from some kind of pathology, but that pathology is different for different density ratios.

For more strongly stratified systems (higher density ratio), 2D simulations produce artificial shear layers that strongly suppress the efficiency of fingering convection, and nonlinearly interact with the latter to produce relaxation oscillations. The emergence of shear at high density ratios in 2D does not come as a surprise, as it was predicted and demonstrated to occur by \citet{Radko2010}. Moreover, a similar phenomenon occurs in 2D Rayleigh-B\'enard convection in vertically-bounded but horizontally-periodic models \citep{Goluskinal14}. The 3D simulations, by contrast, do not show any significant shear. The turbulent fluxes in 3D are significantly higher than in 2D, and are more-or-less steady once the system has reached saturation.

For more weakly stratified systems (lower density ratio), 2D simulations appear to exhibit an inverse energy cascade which results in the progressive coalescence of fingers into larger and larger ones, a phenomenon that is not seen in 3D. As a result, the fluxes at saturation are significantly higher in 2D than in 3D. Furthermore, this inverse cascade seems to prevent the development of the collective instability, which in 3D drives large-scale gravity waves \citep{Medranoal15}.

In short, 2D simulations always fail to reproduce the basic dynamics of fingering convection at low Prandtl number and should therefore never be used in this context. The fact that they appear to be adequate for high-Prandtl number studies \citep{stern2001sfu} is, in the light of our work, somewhat surprising. However, this could be due to the fact that the relative importance of the nonlinear terms in the momentum equation is inversely proportional to the Prandtl number (see equation \ref{eq:nondim}), and that in both cases, the offending artificial behavior discovered in 2D is an inherently nonlinear one. A preliminary study reveals that the transition from dynamically-correct to pathologically-sheared 2D fingering convection occurs for Prandtl numbers around 0.5; in other words, for $\Pr \ge 0.5$, 2D simulations provide qualitatively reasonable results, but for $\Pr < 0.5$, 3D simulations are necessary.

Our findings cast doubt on the 2D fingering flux estimates reported by \citet{denissenkov2010} and \citet{DenissenkovMerryfield2011}. In fact, the simulation snapshots in Figure 4 of \citet{denissenkov2010} are quite similar to our 2D snapshots at $t = 280$ in Figure \ref{fig:2DhiR0snaps}, and clearly show the early stages of the shear development with strongly distorted fingers. Meanwhile, the claim made by \citet{DenissenkovMerryfield2011} that 2D and 3D fluxes at low Prandtl number are very similar could be attributed to the fact that they only integrated their simulations for very short times: we find here that 2D and 3D fluxes are closer to one another before the shear has gained significant amplitude.

We have also found, however, that a fully 3D domain (with equal size in all dimensions) is not necessary to model the correct behavior: a minimum domain width of about 2 wavelengths of the fastest-growing fingering mode is sufficient to eliminate the unwanted behavior, and to yield estimates of the turbulent fluxes that are quantitatively consistent with those obtained in nearly cubic domains. Narrow domains, however, should still be used with caution since they could limit the subsequent development of any large-scale dynamics typically associated with fingering convection (e.g. layer formation and the excitation of large-scale gravity waves).

Interestingly, while there is a fundamental difference between 2D and 3D fingering simulations, \citet{Mollal15} found that this is not the case for simulations in the oscillatory double-diffusive regime (i.e. in the semiconvective regime) at similarly low Prandtl numbers. There, 2D and 3D simulations behave in qualitatively similar ways, and the fluxes extracted in 2D and 3D are within an order of magnitude of each other for all parameters surveyed. This striking difference between the fingering regime and the oscillatory regime raises a number of questions. Why are they so different, given that the only difference in the governing equations is the sign of the background temperature and compositional gradients? Is the ``two-wavelength" guideline we propose for the selection of a minimal 3D domain size in fingering convection exportable to other instabilities? Generally speaking, can one predict ahead of time and without running simulations whether a given set of equations and boundary conditions, solved for a given set of parameters in 2D, will be a good approximation to the full 3D dynamics? These are clearly difficult mathematical questions, but their answers, if they exist, could provide formal guidelines for the general reliability of 2D and narrow-domain 3D simulations.

\acknowledgements

P.G. acknowledges funding from NSF AST-1412951. P. G. and N. B. thank Stephan Stellmach for the development of the code used for this work. All simulations were performed on the UCSC Hyades cluster, purchased thanks to an NSF MRI grant.

\addtocontents{toc}{\protect\vspace*{\baselineskip}}