#  A Bayesian realized threshold measurement GARCH framework for financial
tail risk forecasting

Chao Wang, Richard Gerlach
Discipline of Business Analytics, The University of Sydney

###### Abstract

In this paper, an innovative threshold measurement equation is proposed to be
employed in a Realized-GARCH framework. The proposed framework employs a
nonlinear threshold regression specification to consider the leverage effect
and model the contemporaneous dependence between the observed realized
measures and hidden volatility. A Bayesian Markov Chain Monte Carlo method is
adapted and employed for the model estimation and forecasting, with its
validity assessed via a simulation study. The usefulness of the proposed
measurement equation in a Realized-GARCH model has been evaluated via a
comprehensive empirical study, by forecasting the 1% and 2.5% Value-at-Risk
and Expected Shortfall on six market indices. The proposed framework is shown
to be capable of producing competitive tail risk forecasting results, compared
to the original Realized-GARCH. Especially, the proposed model is favoured
during the high volatility 2008 Global Financial Crisis period.

Keywords: threshold measurement, Realized-GARCH, Markov Chain Monte Carlo,
Value-at-Risk, Expected Shortfall.

##  1 Introduction

A major concern for regulators and owners of financial institutions is the
financial risk management and forecasting. Value-at-Risk (VaR) is one of the
most commonly used risk measures and employed by many financial institutions
as an important risk management tool. VaR represents the market risk as one
number and has become a standard risk measurement metric in recent decades.
Let $\mathcal{I}_{t}$ be the information available at time $t$ and

$F_{t}(r)=Pr(r_{t}\leq r|\mathcal{I}_{t-1})$

be the Cumulative Distribution Function (CDF) of return $r_{t}$ conditional on
$\mathcal{I}_{t-1}$. We assume that $F_{t}(.)$ is strictly increasing and
continuous on the real line $\Re$. Under this assumption, the $\alpha$ level
VaR (quantile) at time $t$ can be defined as:

$Q_{t}=F^{-1}_{t}(\alpha),\qquad 0<\alpha<1.$

However, VaR has been criticized, because it cannot measure the expected loss
for extreme, violating returns and is also not mathematically coherent: i.e.,
it can favour non-diversification. Expected Shortfall (ES), proposed by
Artzner (1997) and Artzner et al. (1999), calculates the expected value of
return being below the quantile of its distribution. Therefore, in recent
years it has become more widely employed for tail risk measurement. The Basel
III Accord, implemented in 2019, places new emphasis on ES . Compared to VaR,
ES possesses a number of attractive properties. For example, ES is a
subadditive risk measure and is mathematically coherent.

Within the same framework as above for defining VaR, the $\alpha$ level ES can
be shown to be equal to the tail conditional expectation of $r_{t}$ (see
Acerbi and Tasche, 2002, among others):

$ES_{t}=E(r_{t}|r_{t}\leq Q_{t,\alpha},\mathcal{I}_{t-1}).$ (1)

It is very important for institutions to have access to highly accurate VaR
and ES forecasts and models, allowing accurate capital allocation, to avoid
both default and over-allocation of funds. As recommended in Basel II and III
Accords, both daily Value-at-Risk (VaR) and Expected Shortfall (ES) are
estimated in this paper, . Basel III Accords recommends the 2.5% level as the
target quantile level which is studied in this paper. A more extreme and
potentially challenging 1% quantile level is also considered.

In parametric VaR and ES estimation and forecasting, producing accurate
volatility estimates and forecasts plays a crucial role. The Autoregressive
Conditional Heteroskedastic (ARCH) and Generalized ARCH (GARCH), proposed in
Engle (1982) and Bollerslev (1986) respectively, gain high popularity in
recent decades.

Black (1976) discovers the now well-known phenomenon of leverage effect
(volatility asymmetry): higher volatility is correlated with negative shocks
in asset returns. Many studies in the literature develop nonlinear GARCH
models to capture it, such as EGARCH (Nelson, 1991) and GJR-GARCH (Glosten et
al., 1993). Poon and Granger (2003) also show that asymmetric volatility
models outperform symmetric models in forecasting asset return volatility.

In addition to EGARCH and GJR-GARCH, another popular type of volatility models
to capture the volatility asymmetry is the threshold autoregressive
specifications, see Tong (1978) and Tong (1990) as examples. A double
threshold GARCH model is developed in Li and Li (1996). Chen et al. (2008)
propose to forecast volatility using threshold heteroskedastic models by
employing the intra-day high-low range.

With the availability of high frequency data, various realized measures have
been proposed and allow potentially more accurate volatility estimation
compared to the daily return, such as Realized Variance (RV) (Andersen and
Bollerslev (1998), Andersen et al. (2003)). Hansen et al. (2012) propose a
Realized-GARCH framework. In addition to the GARCH equation, in Realized-GARCH
a measurement equation that contemporaneously links the unobserved volatility
with a realized measure is included. The Realized-GARCH is shown to be capable
of producing accurate volatility estimation and forecasting results, by
employing RV and other realized measures. The Realized-GARCH model has gained
popularity in the recent decade, with various extensions proposed. A Realized-
EGARCH is developed in Hansen and Huang (2016) which allows for multiple
realized measure to be employed. Huang et al. (2016) develop a long memory
realized heterogeneous autoregressive GARCH (Realized-HAR-GARCH). Chen and
Watanabe (2019) extend the Realized-GARCH model to a threshold framework by
employing a threshold GARCH equation (Realized-Threshold-GARCH). Gerlach and
Wang (2020) propose a semi-parametric realized conditional autoregressive
expectile framework (Realized-CARE).

For all the above mentioned extensions of Realized-GARCH model, in general a
measurement equation that follows the specification in the original Realized-
GARCH model is employed. In this paper, our main contribution is to propose an
innovative threshold measurement equation which can be utilized to capture
volatility asymmetries and the leverage effect. The proposed framework is
motivated by the success of various threshold autoregressive volatility
specifications. The proposed new measurement equation has the same number of
parameters as the original one in the Realized-GARCH and can be employed in
all the above mentioned extended Realized-GARCH models, such as Realized-
EGARCH, Realized-HAR-GARCH, Realized-Threshold-GARCH and Realized-CARE. In our
paper, we mainly evaluate the impact of incorporating this new threshold
measurement equation into the Realized-GARCH, named as Realized-T-M-GARCH. We
would also like to emphasize that the proposed Realized-T-M-GARCH model is
different to the Realized-Threshold-GARCH model in Chen and Watanabe (2019)
which uses a threshold GARCH equation.

Further, an adaptive Bayesian MCMC algorithm is adopted for the estimation of
the proposed model. The validity of the employed MCMC is evaluated via a
simulation study. To evaluate the performance of the proposed Realized-T-M-
GARCH model, employing realized variance as input, the accuracy of the
associated 1% and 2.5% VaR and ES forecasts are assessed via comprehensive
empirical studies. The proposed model produces competitive tail risk
forecasting results compared to the original Realized-GARCH model. Especially,
clear improvements in the risk forecasting accuracy are observed during the
high volatility period, i.e., the 2008 Global Financial Crisis (GFC) period.

The paper is structured as follows: Section 2 reviews the Realized-GARCH model
and proposes the Realized-T-M-GARCH model. The associated likelihood and the
MCMC algorithm for model estimation are presented in Section 3. The simulation
and empirical results are discussed in Section 4 and Section 5 respectively.
Section 6 concludes the paper.

##  2 Model Proposed

###  2.1 Realized-GARCH

The Realized-GARCH model with log specification of Hansen et al. (2012) can be
written as:

$\displaystyle r_{t}$ $\displaystyle=$ $\displaystyle\sqrt{h_{t}}z_{t},$ (2)
$\displaystyle\text{log}(h_{t})$ $\displaystyle=$
$\displaystyle\omega+\beta\text{log}(h_{t-1})+\gamma\text{log}(x_{t-1}),$
$\displaystyle\text{log}(x_{t})$ $\displaystyle=$
$\displaystyle\xi+\varphi\text{log}(h_{t})+\tau_{1}z_{t}+\tau_{2}(z_{t}^{2}-1)+\sigma_{\varepsilon}\varepsilon_{t},$

where $r_{t}=[\text{log}(C_{t})-\text{log}(C_{t-1})]\times 100$ is the
percentage log-return for day $t$, $z_{t}\stackrel{{\scriptstyle\rm
i.i.d.}}{{\sim}}D_{1}(0,1)$, $\varepsilon_{t}\stackrel{{\scriptstyle\rm
i.i.d.}}{{\sim}}D_{2}(0,1)$ and $x_{t}$ is a realized measure, e.g. RV.
$D_{1}(0,1)$ and $D_{2}(0,1)$ indicate distributions that have mean 0 and
variance 1. The three equations in order in model (2) are: the return
equation, the GARCH equation and the measurement equation, respectively. The
measurement equation is an observation equation that captures the dependence
between the latent volatility and the realized measure. The term
$\tau_{1}z_{t}+\tau_{2}(z_{t}^{2}-1)$ is used to capture the leverage effect.

The leverage term in the EGARCH model employs the form
$\tau_{1}z_{t}+\tau_{2}\left(|z_{t}|-E[|z_{t}|]\right)$. Therefore, the
leverage term $\tau_{1}z_{t}+\tau_{2}(z_{t}^{2}-1)$ in Realized-GARCH is a
quadratic variable (constructed from Hermite polynomials) of the version in
EGARCH (see page 270 in Hansen and Huang (2016) for a detailed discussion and
comparison of the two specifications).

Hansen et al. (2012) choose Gaussian errors, e.g. $D_{1}(0,1)=D_{2}(0,1)\equiv
N(0,1)$. Watanabe (2012) allow $D_{1}(0,1)$ to be a standardised Student-t and
skewed Student-t of Fernández and Steel (1998). Student-t is also the choice
of $D_{1}(0,1)$ in Gerlach and Wang (2016). Contino and Gerlach (2017) further
test $D_{2}(0,1)$ as a Student-t distribution, while their findings show that
changing the distribution of $D_{2}(0,1)$ will not significantly affect the
model performance.

###  2.2 Realized threshold measurement GARCH

In a regression model, it is commonly assumed that the coefficients are fixed,
such as the $\text{log}(x_{t})=\xi+\varphi\text{log}(h_{t})$ component in the
measurement equation of the Realized-GARCH. However, in some situations it is
more appropriate to allow the regression coefficients to vary as a function of
time or as a function of some relevant variables. Such regression frameworks
are called switching regression or regime regression model, see Goldfeld and
Quandat (1972), Granger et al. (1993) as examples. The GJR-GARCH utilizes such
framework in modelling the volatility asymmetry and leverage effect, and takes
the following form for the volatility component:

$h_{t}=\omega+\beta h_{t-1}+(\gamma+\alpha I_{t-1})r_{t-1}^{2},\\\ $ (3)

where $I_{t-1}=1$ if $r_{t-1}\leq 0$, otherwise $I_{t-1}=0$. Therefore, the
framework is capable of modelling the volatility asymmetrically according to
the lagged return is positive or negative.

Another popular way of capturing the volatility asymmetry is by incorporating
a threshold GARCH specification (see Li and Li, 1996, among others):

$\displaystyle
h_{t}=\begin{cases}\omega_{1}+\beta_{1}h_{t-1}+\gamma_{1}r_{t-1}^{2},&z_{t-1}\leq
c,\\\
\omega_{2}+\beta_{2}h_{t-1}+\gamma_{2}r_{t-1}^{2},&z_{t-1}>c,\end{cases}$ (4)

where $z_{t}$ is a threshold variable and $c$ is the threshold value.

Motivated by the specifications (3) and (4), we propose an innovative
threshold measurement equation to be employed in the Realized-GARCH framework.
The proposed framework is named as realized threshold measurement GARCH
(Realized-T-M-GARCH) and is shown as:

Realized-T-M-GARCH

$\displaystyle r_{t}$ $\displaystyle=$ $\displaystyle\sqrt{h_{t}}z_{t},$ (5)
$\displaystyle\text{log}(h_{t})$ $\displaystyle=$
$\displaystyle\omega+\beta\text{log}(h_{t-1})+\gamma\text{log}(x_{t-1}),$
$\displaystyle m_{t}$ $\displaystyle=$
$\displaystyle\begin{cases}\xi_{1}+\varphi_{1}\text{log}(h_{t}),&z_{t-1}\leq
c,\\\ \xi_{2}+\varphi_{2}\text{log}(h_{t}),&z_{t-1}>c,\end{cases}$
$\displaystyle\text{log}(x_{t})$ $\displaystyle=$ $\displaystyle
m_{t}+\sigma_{\varepsilon}\varepsilon_{t}.$

In our paper, we choose $z_{t}$ to be self-exciting, that is, $z_{t}=r_{t}$,
and $c=0$, as typical choices in the literature. Therefore, the sign asymmetry
is considered.

Compared to the Realized-GARCH as in (2), the proposed Realized-T-M-GARCH
models the leverage effect in a different (threshold) manner. In addition, the
proposed threshold measurement has exactly the same number of parameters as
the Realized-GARCH.

Although not investigated in our paper, the proposed threshold measurement has
more potential flexibility and could be easily further extended. For example,
the threshold variable $z_{t}$ can be chosen as the realized measure $x_{t}$
and threshold values $c$ can be estimated instead of fixed, thus the size
asymmetry can be also considered. Threshold variable $z_{t}$ can be also
selected as other exogenous economic variables to allow it to be a potentially
more flexible and informative dynamic process.

Further, as mentioned the framework can be easily employed in the other
extensions of the Realized-GARCH framework. Incorporating the threshold
measurement equation into the Realzied-Threshold-GARCH of Chen and Watanabe
(2019), a realized double threshold GARCH can be proposed as below, although
its performance is not investigated in this paper.

Realized-Double-Threshold-GARCH

$\displaystyle r_{t}$ $\displaystyle=$ $\displaystyle\sqrt{h_{t}}z_{t},$ (6)
$\displaystyle\text{log}(h_{t})$ $\displaystyle=$
$\displaystyle\begin{cases}\omega_{1}+\beta_{1}\text{log}(h_{t-1})+\gamma_{1}\text{log}(x_{t-1}),&z_{t-1}\leq
c,\\\
\omega_{2}+\beta_{2}\text{log}(h_{t-1})+\gamma_{2}\text{log}(x_{t-1})&z_{t-1}>c,\end{cases}$
$\displaystyle m_{t}$ $\displaystyle=$
$\displaystyle\begin{cases}\xi_{1}+\varphi_{1}\text{log}(h_{t}),&z_{t-1}\leq
c,\\\ \xi_{2}+\varphi_{2}\text{log}(h_{t}),&z_{t-1}>c,\end{cases}$
$\displaystyle\text{log}(x_{t})$ $\displaystyle=$ $\displaystyle
m_{t}+\sigma_{\varepsilon}\varepsilon_{t},$

Stationarity is an important issue in volatility modelling. As derived in
Hansen et al. (2012) and Gerlach and Wang (2016), by substituting the
threshold measurement equation into volatility equation the required
stationary condition for the Realized-T-M-GARCH model is:

$\beta+\gamma\varphi_{1}<1;\beta+\gamma\varphi_{2}<1.$ (7)

Since the log specification is used, the non-negativity concern of the
volatility is not an issue. Following Watanabe (2012) and Gerlach and Wang
(2016), Student-t distribution is used as the distribution $D_{1}$ of return
equation, while other distributions can be used and tested in the proposed
framework.

##  3 Likelihood and Bayesian Estimation

###  3.1 Likelihood

As in Hansen et al. (2012), when $D_{1}=D_{2}\equiv N(0,1)$ the log-likelihood
function for model (2) is:

$\ell(r,x;\mbox{\boldmath$\theta$})=\underbrace{-\frac{1}{2}\sum_{t=1}^{n}\left[log(2\pi)+log(h_{t})+r_{t}^{2}/h_{t}\right]}_{\ell(r;\mbox{\boldmath$\theta$})}\underbrace{-\frac{1}{2}\sum_{t=1}^{n}\left[log(2\pi)+log(\sigma_{\varepsilon}^{2})+\epsilon_{t}^{2}/\sigma_{\varepsilon}^{2}\right]}_{\ell(x|r;\mbox{\boldmath$\theta$})},$
(8)

where $\epsilon_{t}=log(x_{t})-\xi-\varphi
log(h_{t})-\tau_{1}z_{t}-\tau_{2}(z_{t}^{2}-1)$. The log-likelihood
$\ell(r,x;\mbox{\boldmath$\theta$})$ function equals to the sum of two parts
$\ell(r;\mbox{\boldmath$\theta$})$ and $\ell(x|r;\mbox{\boldmath$\theta$})$
which are derived from the GARCH and measurement equation respectively.

In our paper, we focus on testing the proposed Realized-T-M-GARCH framework
which employs return equation error as $D_{1}\equiv t_{\nu}(0,1)$ and
measurement equation error as $D_{2}\equiv N(0,1)$. $t_{\nu}(0,1)$ represents
a Student-t distribution with $\nu$ degrees of freedom and variance scaled to
1 (by using $\sqrt{\frac{\nu-2}{\nu}}$ factor). The framework is called
Realized-T-M-GARCH-tN.

$\ell(x|r;\mbox{\boldmath$\theta$})$ remains the same as the one in Equation
(8) under the threshold measurement equation specification, as long as the
$D_{2}$ distribution remain unchanged. Therefore, by updating
$\ell(r;\mbox{\boldmath$\theta$})$ according to the Student-t distribution,
the log-likelihood function for the proposed Realized-T-M-GARCH-tN model is:

$\displaystyle\ell(r,x;\theta)$ $\displaystyle=$
$\displaystyle\underbrace{-\sum_{t=1}^{n}\left[A(\nu)+\frac{1}{2}\log(h_{t})+\frac{\nu+1}{2}\left(1+\frac{r_{t}^{2}}{h_{t}(\nu-2)}\right)\right]}_{\ell(r;\theta)}$
$\displaystyle\underbrace{-\frac{1}{2}\sum_{t=1}^{n}\left[log(2\pi)+log(\sigma_{\varepsilon}^{2})+\epsilon_{t}^{2}/\sigma_{\varepsilon}^{2}\right]}_{\ell(r|x;\theta)}$

where $\epsilon_{t}=x_{t}-\xi_{1}-\varphi_{1}h_{t}$ when $r_{t-1}\leq 0$,
$\epsilon_{t}=x_{t}-\xi_{2}-\varphi_{2}h_{t}$ when $r_{t-1}>0$, and
$A(\nu)=-\log\left(\Gamma\left(\frac{\nu+1}{2}\right)\right)+\frac{1}{2}\log(\pi(\nu-2))+\log\left(\Gamma\left(\frac{\nu}{2}\right)\right)$.
The parameter vector to be estimated is now
$\mbox{\boldmath$\theta$}=(\omega,\beta,\gamma,\xi_{1},\varphi_{1},\xi_{2},\varphi_{2},\sigma_{\varepsilon},\nu)^{{}^{\prime}}$,
under the constraints in Equation (7). $\nu>4$ is further restricted to ensure
the first four moments of the error distribution are finite.

For comparison purpose, we also incorporate Student-t distribution as $D_{1}$
and Normal distribution as $D_{2}$ in the Realized-GARCH (Realized-GARCH-tN).
It has identical log-likelihood as the Realized-T-M-GARCH-tN framework.

###  3.2 Bayesian Estimation

Motivated by the MCMC results in Gerlach and Wang (2016) and Gerlach and Wang
(2020), an adaptive Markov Chain Monte Carlo (MCMC) procedure is employed for
the estimation and forecasting of both Realized-T-M-GARCH-tN and ealized-
GARCH-tN, also to make the later model performance comparison a fair one. The
likelihood for both models, as in Equation (3.1), involves 9 unknown
parameters.

The motivations of employing the Bayesian MCMC approach (instead of maximum
likelihood estimation (MLE)) are discussed in Gerlach and Wang (2016). For
example, the estimation of Realized-T-M-GARCH and Realized-GARCH require
constrained MLE, to ensure stationarity, and this may cause issues in the
optimisation and the standard error calculation (Silvapulle and Sen, 2005).
Therefore, in this paper we have adapted and extended the adaptive MCMC method
in Gerlach and Wang (2020) which contains two steps: a burn-in step and a
“independent” Metropolis-Hastings (IMH) step.

First, the parameter blocking is employed for the MCMC estimation of Realized-
T-M-GARCH-tN. Three blocks are chosen as:
$\mbox{\boldmath$\theta_{1}$}=(\omega,\beta,\gamma,\varphi_{1},\varphi_{2})^{{}^{\prime}}$;
$\mbox{\boldmath$\theta_{2}$}=(\xi_{1},\xi_{2},\sigma_{\varepsilon})^{{}^{\prime}}$
and $\mbox{\boldmath$\theta_{3}$}=(\nu)$. The choice is motivated by the face
that parameters within the same block are more strongly correlated, in the
posterior (or likelihood), than those between blocks. For example, the
stationarity condition causes correlation between iterates of
$\beta,\gamma,\varphi_{1},\varphi_{2}$, thus they are kept together in one
block. Similarly, three blocks for the Realized-GARCH-tN model are also used:
$\mbox{\boldmath$\theta_{1}$}=(\omega,\beta,\gamma,\varphi)^{{}^{\prime}}$;
$\mbox{\boldmath$\theta_{2}$}=(\xi,\tau_{1},\tau_{2},\sigma_{\varepsilon})^{{}^{\prime}}$
and $\mbox{\boldmath$\theta_{3}$}=(\nu)$. Uninformative prior are chosen for
both models over the possible stationarity region, i.e.,
$\pi(\mbox{\boldmath$\theta$})\propto I(A)$, which is a flat prior for
$\theta$ over the region $A$ which satisfies the stationarity condition in
Equation (7).

For the burn-in period, a Metropolis algorithm (Metropolis _et al._ , 1953)
employing a mixture of 3 Gaussian proposal distributions, with a random walk
mean vector, is employed for each block of parameters. In addition, an
iterative “epoch” method in Chen et al. (2017) is employed, aim to enhance the
convergence of MCMC chains. The proposal variance-covariance (var-cov) matrix
of each block in each mixture element is set as $C_{i}\Sigma_{1}$, where
$C_{1}=1;C_{2}=100;C_{3}=0.01$, with $\Sigma_{1}$ initially set to
$\frac{2.38}{\sqrt{(d_{i})}}I_{d_{i}}$, where $d_{i}$ is the dimension of the
block ($i$) of parameters being generated, and $I_{d_{i}}$ is the identity
matrix of dimension $d_{i}$. This proposal var-cov matrix is subsequently
tuned, aiming towards a target acceptance rate of $23.4\%$ (if $d_{i}>4$, or
$35\%$ if $2\leq d_{i}\leq 4$, or $44\%$ if $d_{i}=1$), as standard, via the
algorithm of Gelman et al. (1997). After running the first epoch with 20,000
iterations, the var-cov matrix of each parameter block is calculated after
discarding the first 2,000 iterates. Such updated var-cov matrix is then
employed in the proposal distribution of next epoch. After running each epoch,
the standard deviations of each parameter chain are also calculated and
compared to that of the previous epoch. Such epoch process is repeated until
the mean absolute percentage change of these standard deviations is less than
10%, which takes approximately three to four epochs in both the simulation and
empirical study.

In the IMH step (10,000 iterations), again a mixture of three Gaussian
proposal distributions is utilised for each parameter block. The sample mean
vector of the last epoch iterates (after discarding the first 2,000 iterates)
in the burn-in period is used as the mean vector of the IMH step. For var-cov
matrix, after discarding the first 2,000 iterates, the sample covariance
matrix of the last epoch iterates for each block is calculated as
$\Sigma_{2}$. Then the proposal variance-covariance matrix in each element is
calculated as $C_{i}\Sigma_{2}$, where $C_{1}=1;C_{2}=100;C_{3}=0.01$.

Lastly, all the IMH iterates (still discarding the first 2,000 iterates) are
employed to calculate the VaR and ES forecasts, the posterior mean of which is
used as the final tail risk forecasts.

##  4 Simulation Study

A simulation study is designed to illustrate the validity of the adapted MCMC,
in terms of parameter estimation, VaR and ES forecasting accuracy.

1000 replicated data sets of size $n=1900$ (chosen based on empirical study
in-sample size) are simulated from the following simulation model which
incorporates the Realized-T-M-GARCH-tN specification.

Simulation Model

$\displaystyle r_{t}$ $\displaystyle=$
$\displaystyle\sqrt{h_{t}}z_{t},\,\,z_{t}\sim t_{\nu=10}(0,1)$ (10)
$\displaystyle\text{log}(h_{t})$ $\displaystyle=$ $\displaystyle
0.1+0.65\text{log}(h_{t-1})+0.3\text{log}(x_{t-1}),$ $\displaystyle m_{t}$
$\displaystyle=$
$\displaystyle\begin{cases}-0.2+0.95\text{log}(h_{t}),&r_{t-1}\leq 0,\\\
-0.5+0.92\text{log}(h_{t}),&r_{t-1}>0,\end{cases}$
$\displaystyle\text{log}(x_{t})$ $\displaystyle=$ $\displaystyle
m_{t}+0.6\varepsilon_{t}\,,$

The true values of the parameters in the simulation model are selected based
on the parameter estimates of the empirical study. The estimated parameter
values of Realized-T-M-GARCH-tN are in general consistent with the ones in the
Realized-GARCH, while distinctive behaviours are observed for the parameters
in the measurement equation.

Specifically, in a Realized-GARCH the $\gamma$ parameter is typically
estimated to be between 0.3 and 0.55 in the empirical study. This parameter
may be compared with $\alpha$ in a conventional GARCH model, which measures
the coefficient associated with volatility estimator (square return in GARCH).
However, with a more efficient and informative realized measure employed, the
estimated $\gamma$ parameter is in general greater than estimated $\alpha$ in
GARCH. Under the Realized-T-M-GARCH framework, we have similar observations.
Thus, the true values of $\gamma$ in the simulation model is selected as 0.3.

In a Realized-GARCH, the estimates of $\varphi$ are close to unity, which
suggests that the realized measure $x_{t}$ is roughly proportional to the
conditional variance of daily returns. The parameter $\xi$ in a Realized-GARCH
is always negative. This suggests that the incorporated realized measures are
roughly proportional to the conditional variance, with a small negative
correction (given by estimates of $\xi$). One potential reason for this is
that the returns employed in this paper are close-to-close and include
overnight price movements, but the realized variance is only measured when the
market is open which may slightly underestimate true volatility on average
(downside bias). The leverage effect is prominent in stock indices, with the
estimates of $\tau_{1}$ always as negative in a Realized-GARCH.

To support the discussion on the threshold measurement equation true parameter
value selection of the Realized-T-M-GARCH-tG simulation model, we present the
S&P500 out-of-sample $\xi_{1}$&$\xi_{2}$, $\varphi_{1}$&$\varphi_{2}$
parameter estimates in Figure 1. The details of empirical results will be
discussed later in the empirical section. In general, Figure 1 shows that the
$\xi_{1}$&$\xi_{2}$, $\varphi_{1}$&$\varphi_{2}$ in the threshold measurement
display distinctive behaviours with different estimated values.

Figure 1:  S&P 500 out-of-sample threshold measurement equation parameter
estimates and persistent level.

More specifically, as in the second plot in Figure 1 the $\xi_{2}$ estimated
value is tend to be more negative (when lagged return is positive,
$r_{t-1}>0$) than the $\xi_{1}$ estimated value ($r_{t-1}\leq 0$). Such
results lend evidence on the fact of leverage effect which is successfully
captured by the $\xi_{1}$ and $\xi_{2}$ parameters, with more negative bias
correction is required when the lagged return is positive. Therefore, in the
simulation model the true values for $\xi_{1}$ and $\xi_{2}$ are selected as
-0.2 and -0.5 respectively. This finding is actually consistent with the
Realized-GARCH $\tau_{1}$ estimate which is always negative, so when a return
(or correspondingly the error term $z_{t}$) is positive there is more negative
bias correction required.

In addition, comparing to the Realized-GARCH which has one regression
coefficient $\varphi$ that links the observed realized measure $x_{t}$ and
hidden conditional variance $h_{t}$. The proposed Realized-T-M-GARCH includes
two separate coefficients $\varphi_{1}$ and $\varphi_{2}$ in two regimes which
can potentially model the relationship between $x_{t}$ and $h_{t}$ more
flexibly. The estimates of $\varphi_{1}$ and $\varphi_{2}$ in the empirical
study provide potential support on this. Both parameter are estimated to be
close to 1, while different values are observed at different forecasting
steps, as in the third plot in Figure 1. The true parameter values for
$\varphi_{1}$ and $\varphi_{2}$ are chosen as 0.95 and 0.92 respectively.
Lastly, although the persistence levels of two regimes (fourth plot in Figure
1), $\beta+\gamma\varphi_{1}$ and $\beta+\gamma\varphi_{2}$, are both very
close to 1, distinctive behaviours are again observed during the out-of-sample
period.

Overall, the true parameter values in the simulation model represent typical
results observed in the empirical study. In the simulation model, $r_{t}$ is
analogous to a daily return and $x_{t}$ is analogous to the daily realized
measure.

Simulation results of both 1% and 2.5% quantile levels are presented. All
initial parameter values are arbitrarily set equal to $0.25$ to start the MCMC
chains.

The “True” one-step-ahead level $\alpha$ VaR forecast from the above
simulation model is calculated as:

$Q_{t+1}=\sqrt{h_{t+1}}t_{\nu}^{-1}(\alpha)\sqrt{\frac{\nu-2}{\nu}},$

where $t_{\nu}^{-1}(\alpha)$ is the inverse of Student-t’s CDF with the $\nu$
degrees of freedom on quantile level $\alpha$. $\sqrt{\frac{\nu-2}{\nu}}$ is
used to scale the Student-t distribution to variance 1. ES forecast from the
same model is calculated as:

$\text{ES}_{t+1}=-\sqrt{h_{t+1}}\left(\frac{g_{\nu}(t_{\nu}^{-1}(\alpha))}{\alpha}\right)\left(\frac{\nu+(t_{\nu}^{-1}(\alpha))^{2}}{\nu-1}\right)\sqrt{\frac{\nu-2}{\nu}},$

where $g_{\nu}$ is the Student-t Probability Density Function (PDF).

These “True” VaR and ES forecasts are calculated for each data set. The
averages of the these “True” forecasts, over the 1000 data sets, are given in
the “True” column of Table 1.

The simulation results are summarised in Table 1. The Mean column shows the
average parameter estimates and average VaR&ES forecasts across all 1000
simulated data sets. The Root Mean Squared Error (RMSE) between the parameter
estimates and parameter true values (and VaR and ES forecasts and their “True”
values) are also shown. In general, the MCMC algorithm produces accurate
parameter estimates and tail-risk forecasts in terms of bias (difference
between True and Mean columns) and precision (RMSE). Specifically, both 1% and
2.5% VaR&ES forecasts are close to their true values with bias less than 0.05
and RMSE less than 0.15. This means the employed MCMC procedure could produce
accurate parameter estimates and tail risk forecasts.

Table 1:  Summary statistics for the MCMC estimator of the Realized-T-M-GARCH model, data simulated from Model 1. $n=1900$ |  | MCMC
---|---|---
Parameter | True | Mean | RMSE
$\omega$ | 0.1000 | 0.0993 | 0.0169
$\beta$ | 0.6500 | 0.6426 | 0.0236
$\gamma$ | 0.3000 | 0.3021 | 0.0259
$\xi_{1}$ | -0.2000 | -0.1983 | 0.0420
$\varphi_{1}$ | 0.9500 | 0.9635 | 0.0811
$\xi_{2}$ | -0.5000 | -0.4935 | 0.0426
$\varphi_{2}$ | 0.9200 | 0.9320 | 0.0835
$\sigma_{\varepsilon}$ | 0.6000 | 0.6011 | 0.0100
$\nu$ | 10.0000 | 12.3515 | 3.9844
$\text{1\% VaR}_{t+1}$ | -2.4576 | -2.4424 | 0.0905
$\text{2.5\% VaR}_{t+1}$ | -1.9813 | -1.9762 | 0.0616
$\text{1\% ES}_{t+1}$ | -2.9907 | -2.9589 | 0.1456
$\text{2.5\% VaR}_{t+1}$ | -2.5068 | -2.4886 | 0.0996

Further, since the main development of the paper is the threshold measurement
equation, Figure 2 presents the histograms of parameter estimates of all 1000
simulated data sets for $\xi_{1}$&$\xi_{2}$ and $\varphi_{1}$&$\varphi_{2}$.
As can be seen, although parameters $\xi_{1}$&$\xi_{2}$ and
$\varphi_{1}$&$\varphi_{2}$ have different true values, the adapted MCMC
algorithm is capable of estimating all four parameters accurately, with the
RMSE values around 0.04 for $\xi_{1}$&$\xi_{2}$, and around 0.08 for
$\varphi_{1}$&$\varphi_{2}$ respectively.

Figure 2:  Histograms of 1000 parameter estimates for $\xi_{1}$&$\xi_{2}$ and
$\varphi_{1}$&$\varphi_{2}$ in the threshold measurement equation. True
vertical line represents parameter true value from the simulation model. Mean
vertical line represents the average of 1000 parameter estimates from MCMC.

##  5 Empirical Study

###  5.1 Data and forecasting setup

Six market indices are assessed in the empirical section, including S&P500,
NASDAQ (both US), FTSE 100 (UK), DAX (Germany), SMI (Swiss) and ASX200
(Australia), with time period from January 2000 to December 2015. The high
frequency closing prices, observed at 5-minute intervals within trading hours,
are downloaded from Thomson Reuters Tick History. The 5-minute data are
employed to calculate the daily RV. The daily closing prices are also
collected and used to calculate the daily return.

As discussed in Section 4, the daily returns employed are close-to-close,
including the overnight price movement, while the RV is only calculated while
the market is open. Therefore, the employed realized measures may have some
downside bias for the true volatility of close-to-close returns, which can be
captured and corrected by the $\xi_{1}$ and $\xi_{2}$ parameters in the
Realized-T-M-GARCH framework.

A fixed size in-sample data is employed for estimation, combined with a
rolling window approach, to produce each one-step-ahead VaR and ES forecast.
Table 2 reports the in-sample size for each series, which differs due to
different non-trading days occurring in each market.

Two forecasting studies with different out-of-sample sizes are conducted. The
first study aims to assess the performance of the models for the 2008 GFC high
volatility period, thus the initial date of the out-of-sample forecasting
period is chosen as January 2008. Then for each index the out-of-sample size
$m$ is chosen as 400, meaning that the end of the forecasting period is around
August 2009.

An eight year out-of-sample period is employed in the second forecasting study
with the start date of the out-of-sample still chosen as Jan 2008 and out-of-
sample size $m$ as 2000. Therefore, the end of the forecasting period is
around December 2015.

Both daily one-step-ahead VaR and ES forecasts are considered for the returns
on the six indices, using $\alpha=1\%;2.5\%$.

The proposed Realized-T-M-GARCH-tG and Realized-GARCH-tG, both estimated with
MCMC, are included in the forecasting study to compare their performance.

In addition, for comparison purpose EGARCH and GJR-GARCH, both with Student-t
distribution, are also included. Further, a GARCH with historical simulation
(GARCH-HS) approach is also included, where either a GJR-GARCH-t or a EGARCH-t
is fit to the in-sample data, then standardised VaR and ES are estimated via
historical simulation from the sample of returns (e.g. $r_{1},\ldots,r_{n}$),
standardised by dividing their GARCH-estimated conditional standard deviation,
i.e., $r_{t}/\sqrt{\hat{h}_{t}}$. Then, forecasts of VaR, ES are found by
multiplying the standardised VaR and ES estimates by the forecast
$\sqrt{\hat{h}_{n+1}}$ from the GJR-GARCH-t or the EGARCH-t model. All these
GARCH type models are estimated by MLE.

###  5.2 VaR forecasting

Since quantiles are elicitable, as defined in Gneiting (2011), and the
standard quantile loss function is strictly consistent, which means the
expected quantile loss is a minimum at the true quantile series. In this
section, the quantile loss over the out-of-sample period is used to compare
the VaR forecast accuracy of the competing models. The most accurate VaR
forecasting model is expected to produce the minimized aggregated quantile
loss function, given as in Equation (11):

$\sum_{t=n+1}^{n+m}(\alpha-I(r_{t}<Q_{t}))(r_{t}-Q_{t})\,\,,$ (11)

where $n$ is the in-sample size and $m$ is the out-of-sample size.
$\widehat{Q}_{n+1},\ldots,\widehat{Q}_{n+m}$ is a series of quantile forecasts
at levels $\alpha=1\%;2.5\%$ for the returns $r_{n+1},\ldots,r_{n+m}$.

Values of the out-of-sample quantile loss are presented in Table 2 and 3. The
average loss is included in the Avg Loss column. The average rank based on
ranks of quantile loss across six markets is calculated and shown in the Avg
Rank column. Box indicates the favoured model and dashed box indicates the 2nd
ranked model based on the average loss and rank.

For the GFC forecasting study, Table 2 shows that the proposed Realized-T-M-
GARCH-tG framework is characterized by very competitive performances, on both
1% and 2.5% quantile levels. On the 1% level, Realized-T-M-GARCH-tG ranks as
the best on average (1.33), with smallest average quantile loss produced
(25.2). Its performance is followed by Realized-GARCH-tG which has the same
average rank as EGARCH-t-HS. On the 2.5% level, the best ranked model (2.00)
with the smallest average quantile loss (52.4) is still Realized-T-M-GARCH-tG.
The model with the second smallest loss is Realized-GARCH-tG, while the second
best ranked model is EGARCH-t-HS. The performance of EGARCH-t and GJR-GARCH-t
is in general close. Lastly, the effectiveness of incorporating realized
measure in the quantile forecasting is clear, with the Realized-T-M-GARCH-tG
and Realized-GARCH-tG in general better ranked with smaller loss produced.

Table 2:  For the GFC study, quantile loss function values across the markets.
Out-of-sample size $m=400$. $\alpha=1\%;2.5\%$ .

Model | S&P500 | NASDAQ | FTSE | DAX | SMI | ASX200 | Avg Loss | Avg Rank
---|---|---|---|---|---|---|---|---
$\alpha=1\%$ |  |  |  |  |  |  |  |
EGARCH-t | 27.3 | 32.8 | 26.1 | 26.3 | 26.9 | 21.5 | 26.8 | 4.33
GJR-GARCH-t | 26.0 | 32.0 | 26.8 | 27.1 | 27.8 | 22.3 | 27.0 | 4.83
EGARCH-t-HS | 26.1 | 32.4 | 24.8 | 26.5 | 26.1 | 20.7 | 26.1 |  3.00
GJR-GARCH-t-HS | 26.2 | 32.2 | 25.7 | 27.0 | 27.0 | 22.0 | 26.7 | 4.50
Realized-GARCH-tG | 26.7 | 30.0 | 24.4 | 26.8 | 26.3 | 20.0 |  25.7 |  3.00
Realized-T-M-GARCH-tG | 25.7 | 29.2 | 24.3 | 26.6 | 25.3 | 20.0 | 25.2 | 1.33
$\alpha=2.5\%$ |  |  |  |  |  |  |  |
EGARCH-t | 59.8 | 62.8 | 54.7 | 53.4 | 51.4 | 47.8 | 55.0 | 4.50
GJR-GARCH-t | 55.4 | 62.0 | 54.9 | 54.0 | 52.3 | 47.2 | 54.3 | 4.50
EGARCH-t-HS | 57.1 | 62.7 | 52.3 | 52.6 | 50.8 | 45.5 | 53.5 |  3.00
GJR-GARCH-t-HS | 55.1 | 62.0 | 53.1 | 53.8 | 52.2 | 46.1 | 53.7 | 3.33
Realized-GARCH-tG | 56.2 | 59.4 | 51.2 | 55.3 | 52.6 | 42.6 |  52.9 | 3.67
Realized-T-M-GARCH-tG | 54.9 | 58.7 | 51.1 | 55.3 | 51.7 | 42.5 | 52.4 | 2.00
Out-of-sample $m$ | 400 | 400 | 400 | 400 | 400 | 400 |  |
In-sample $n$ | 2018 | 2003 | 2081 | 2066 | 2033 | 1986 |  |

_Note_ :Based on average loss and rank, box indicates the favoured model,
dashed box indicates the 2nd ranked model.

Table 3 includes the quantile forecasting results for the eight year out-of-
sample study. For the more extreme 1% quantile level, the first and second
best performed models are still Realized-T-M-GARCH-tG and Realized-GARCH-tG.
However, for the 2.5% level, the best performing model is EGARCH-t-HS with
average loss 162.4. Its performance is closely followed by Realized-T-M-GARCH-
tG with average loss 162.5. The GJR-GARCH-t-HS model also produces competitive
results and is second best ranked (together with Realized-T-M-GARCH-tG).

To conclude, for both forecasting studies the proposed Realized-T-M-GARCH-tG
model is characterised by very competitive quantile loss results. This
demonstrates the validity and effectiveness of incorporating the threshold
measurement equation to consider the leverage effect and forecast VaR
accurately.

Table 3:  Quantile loss function values across the markets. Out-of-sample size
$m=2000$. $\alpha=1\%;2.5\%$.

Model | S&P500 | NASDAQ | FTSE | DAX | SMI | ASX200 | Avg Loss | Avg Rank
---|---|---|---|---|---|---|---|---
$\alpha=1\%$ |  |  |  |  |  |  |  |
EGARCH-t | 77.4 | 88.8 | 72.5 | 86.3 | 79.9 | 64.3 | 78.2 | 4.33
GJR-GARCH-t | 74.9 | 86.6 | 73.3 | 87.8 | 82.6 | 64.8 | 78.3 | 5.00
EGARCH-t-HS | 76.7 | 88.3 | 71.3 | 86.9 | 78.9 | 64.0 | 77.7 | 3.17
GJR-GARCH-t-HS | 75.2 | 86.5 | 71.9 | 88.0 | 81.7 | 64.6 | 78.0 | 4.17
Realized-GARCH-tG | 74.9 | 84.5 | 72.2 | 84.3 | 80.8 | 62.9 |  76.6 |  2.50
Realized-T-M-GARCH-tG | 74.7 | 84.1 | 72.2 | 84.3 | 80.5 | 62.8 | 76.4 | 1.83
$\alpha=2.5\%$ |  |  |  |  |  |  |  |
EGARCH-t | 166.9 | 184.0 | 155.2 | 181.6 | 159.3 | 140.3 | 164.6 | 4.50
GJR-GARCH-t | 162.9 | 183.1 | 156.4 | 182.9 | 159.4 | 141.7 | 164.4 | 4.83
EGARCH-t-HS | 163.7 | 182.3 | 153.0 | 179.5 | 157.5 | 138.4 | 162.4 | 2.50
GJR-GARCH-t-HS | 161.0 | 180.6 | 154.5 | 180.7 | 159.3 | 139.8 | 162.6 |  2.83
Realized-GARCH-tG | 162.0 | 178.7 | 154.1 | 185.1 | 161.8 | 134.9 | 162.8 | 3.50
Realized-T-M-GARCH-tG | 161.1 | 178.6 | 154.0 | 185.2 | 161.2 | 134.9 |  162.5 |  2.83
Out-of-sample $m$ | 2000 | 2000 | 2000 | 2000 | 2000 | 2000 |  |
In-sample $n$ | 2018 | 2003 | 2081 | 2066 | 2033 | 1986 |  |

_Note_ :Based on average loss and rank, box indicates the favoured model,
dashed box indicates the 2nd ranked model.

###  5.3 ES forecasting

The same set of models is employed to generate one-step-ahead forecasts of 1%
and 2.5% ES during the forecast period for all six series.

To evaluate the proposed Realized-T-M-GARCH, we assess the ability of the
different models under comparison to forecast VaR and ES jointly, employing a
strictly consistent VaR and ES joint loss function.

Fissler and Ziegel (2016) find the class of jointly consistent scoring
functions for VaR and ES i.e., their expectations are uniquely minimized by
the true VaR and ES series. The general form of this functional family is:

$\displaystyle S_{t}(r_{t},Q_{t},ES_{t})$ $\displaystyle=$
$\displaystyle(I_{t}-\alpha)G_{1}(Q_{t})-I_{t}G_{1}(r_{t})+G_{2}(ES_{t})\left(ES_{t}-Q_{t}+\frac{I_{t}}{\alpha}(Q_{t}-r_{t})\right)$
(12) $\displaystyle-$ $\displaystyle H(ES_{t})+a(r_{t})\,,$

where $G_{1}(.)$ is increasing, $G_{2}(.)$ is strictly increasing and strictly
convex, $G_{2}=H^{{}^{\prime}}$ and $\lim_{x\to-\infty}G_{2}(x)=0$ and
$a(\cdot)$ is a real-valued integrable function.

As presented in Taylor (2019), assuming $r_{t}$ to have zero mean, making the
choices: $G_{1}(x)=0$, $G_{2}(x)=-1/x$, $H(x)=-\text{log}(-x)$ and
$a=1-\text{log}(1-\alpha)$, which satisfy the required criteria, returns the
scoring function:

$\displaystyle
S_{t}(r_{t},Q_{t},ES_{t})=-\text{log}\left(\frac{\alpha-1}{\text{ES}_{t}}\right)-{\frac{(r_{t}-Q_{t})(\alpha-I(r_{t}<Q_{t}))}{\alpha\text{ES}_{t}}}.$
(13)

Taylor (2019) refers to Equation (13) as the Asymmetric Laplace (AL) log score
which is used to jointly assess the VaR and ES forecasting accuracy.

First, Figure 3 shows the S&P500 1% ES forecasts from EGARCH-t, Realized-
GARCH-tG and Realized-T-M-GARCH-tG during the GFC period. As can be seen, the
ES forecasts from Realized-T-M-GARCH-tG and Realized-GARCH-tG present some
distinctive behaviours during the highly volatile period, i.e., around Oct
2008. To make a more in-depth comparison of these models, Figure 4 presents
the S&P 500 1% VaR&ES AL joint loss (log-score) values for each time step
across the out-of-sample period of the GFC study. In general, the Realized-T-
M-GARCH-tG VaR&ES forecasts are characterized by smaller joint loss values
than the ones from the Realized-GARCH-tG, such as Oct 2008 and Dec 2008. In
addition, there are two joint loss value jumps for the Realized-GARCH-tG
during the first half of 2009, while such jumps are not observed for the
Realized-T-M-GARCH-tG model. Such reduced VaR&ES joint loss values, especially
during the high volatility period, reflect the potential additional tail risk
forecasting efficiency that can be gained from the threshold measurement
specification in a Realized-GARCH framework. These observations are also
consistently presented across different indices.

Figure 3:  S&P 500 1% ES forecasts from EGARCH-t, Realized-GARCH-tG and
Realized-T-M-GARCH-tG during the GFC period. Figure 4:  S&P 500 1% VaR&ES
forecasts AL joint loss values from EGARCH-t, Realized-GARCH-tG and Realized-
T-M-GARCH-tG during the GFC period.

As a more comprehensive comparison, Table 4 and Table 5 report, for each model
and index, the value of the loss function in Equation (13) aggregated over the
out-of-sample period: $\mathbf{S}=\sum_{t=n+1}^{n+m}S_{t}$, with $m=400$ and
$m=2000$ respectively.

Regarding the joint loss of the GFC study, Realized-T-M-GARCH-tG is again the
best ranked model with the smallest average loss produced, on both 1% and 2.5%
quantile levels. The second best performing model is Realized-GARCH-t.
Comparing models with and without high frequency information, the performance
of Realized-T-M-GARCH-tG and Realized-GARCH-tG is clearly favoured.

Table 4:  For the GFC study, VaR and ES joint loss function values across the
markets. Out-of-sample size $m=400$. $\alpha=1\%;2.5\%$.

Model | S&P500 | NASDAQ | FTSE | DAX | SMI | ASX200 | Avg Loss | Avg Rank
---|---|---|---|---|---|---|---|---
$\alpha=1\%$ |  |  |  |  |  |  |  |
EGARCH-t | 1167.6 | 1244.6 | 1205.5 | 1161.1 | 1135.5 | 1066.2 | 1163.4 | 4.50
GJR-GARCH-t | 1123.9 | 1228.1 | 1208.2 | 1186.8 | 1170.6 | 1088.8 | 1167.7 | 5.50
EGARCH-t-HS | 1134.0 | 1226.3 | 1160.9 | 1156.3 | 1106.0 | 1039.6 | 1137.2 | 3.00
GJR-GARCH-t-HS | 1119.8 | 1226.4 | 1164.1 | 1174.5 | 1139.1 | 1068.3 | 1148.7 | 3.83
Realized-GARCH-tG | 1120.0 | 1182.7 | 1103.3 | 1176.2 | 1104.2 | 1028.0 |  1119.1 |  2.50
Realized-T-M-GARCH-tG | 1102.0 | 1173.4 | 1101.1 | 1179.8 | 1083.7 | 1028.0 | 1111.3 | 1.67
$\alpha=2.5\%$ |  |  |  |  |  |  |  |
EGARCH-t | 1114.2 | 1131.5 | 1108.8 | 1068.2 | 1031.9 | 1030.3 | 1080.8 | 4.83
GJR-GARCH-t | 1078.6 | 1124.4 | 1107.0 | 1082.6 | 1047.7 | 1027.9 | 1078.1 | 4.67
EGARCH-t-HS | 1085.1 | 1126.6 | 1073.3 | 1058.4 | 1019.4 | 995.1 | 1059.7 | 3.00
GJR-GARCH-t-HS | 1069.3 | 1122.5 | 1073.9 | 1076.9 | 1037.7 | 1007.2 | 1064.6 | 3.67
Realized-GARCH-tG | 1067.6 | 1091.4 | 1038.6 | 1088.3 | 1035.4 | 975.6 |  1049.5 |  2.83
Realized-T-M-GARCH-tG | 1056.7 | 1085.7 | 1037.4 | 1091.1 | 1019.5 | 975.1 | 1044.3 | 2.00
Out-of-sample $m$ | 400 | 400 | 400 | 400 | 400 | 400 |  |
In-sample $n$ | 2018 | 2003 | 2081 | 2066 | 2033 | 1986 |  |

_Note_ :Based on average loss and rank, box indicates the favoured model,
dashed box indicates the 2nd ranked model.

With respect to the joint loss values produced with longer forecasting
horizon, on the 1% quantile level the preferred model with the smallest
average loss and best rank is Realized-GARCH-tG, which is closely followed by
the Realized-T-M-GARCH-tG framework. With $\alpha=2.5\%$ which is a less
extreme quantile level than 1%, Realized-T-M-GARCH-tG generates competitive
results compared to Realized-GARCH-tG, with both models ranked as 2.67. The
joint loss values of two models are very close as well, 4244.0 and 4242.5
respectively. The best ranked model is GJR-GARCH-t-HS (2.00).

Overall, the joint loss results lend further support on the usefulness of the
proposed Realized-T-M-GARCH framework. The model produces competitive results,
comparing to the competing models, for both forecasting studies and both
quantile levels. Especially, for the highly volatile GFC period, the Realized-
T-M-GARCH framework is clearly preferred, with potential reasons as discussed
when presenting Figure 3 and 4.

Table 5:  VaR and ES joint loss function values across the markets. Out-of-
sample size $m=2000$. $\alpha=1\%;2.5\%$.

Model | S&P500 | NASDAQ | FTSE | DAX | SMI | ASX200 | Avg Loss | Avg Rank
---|---|---|---|---|---|---|---|---
$\alpha=1\%$ |  |  |  |  |  |  |  |
EGARCH-t | 4587.6 | 4837.2 | 4537.7 | 4949.5 | 4689.9 | 4277.6 | 4646.6 | 5.67
GJR-GARCH-t | 4459.0 | 4743.6 | 4486.0 | 4981.8 | 4778.7 | 4245.0 | 4615.7 | 4.83
EGARCH-t-HS | 4539.0 | 4782.7 | 4479.9 | 4924.8 | 4604.2 | 4258.8 | 4598.2 | 3.83
GJR-GARCH-t-HS | 4436.8 | 4717.4 | 4422.8 | 4947.4 | 4684.5 | 4225.1 | 4572.3 | 3.33
Realized-GARCH-tG | 4415.8 | 4715.0 | 4420.3 | 4822.6 | 4630.5 | 4149.5 | 4525.6 | 1.50
Realized-T-M-GARCH-tG | 4433.0 | 4706.7 | 4421.0 | 4831.6 | 4646.0 | 4148.9 |  4531.2 |  1.83
$\alpha=2.5\%$ |  |  |  |  |  |  |  |
EGARCH-t | 4290.6 | 4517.6 | 4192.1 | 4585.2 | 4245.2 | 3990.4 | 4303.5 | 5.33
GJR-GARCH-t | 4239.2 | 4490.2 | 4183.8 | 4601.4 | 4253.6 | 4009.1 | 4296.2 | 5.50
EGARCH-t-HS | 4225.1 | 4464.5 | 4146.9 | 4542.8 | 4185.5 | 3953.9 | 4253.1 | 2.83
GJR-GARCH-t-HS | 4174.1 | 4428.8 | 4131.3 | 4556.6 | 4216.0 | 3962.5 | 4244.9 | 2.00
Realized-GARCH-tG | 4183.9 | 4429.4 | 4144.5 | 4582.9 | 4226.2 | 3888.0 | 4242.5 |  2.67
Realized-T-M-GARCH-tG | 4188.0 | 4426.3 | 4144.3 | 4591.5 | 4225.6 | 3888.3 |  4244.0 |  2.67
Out-of-sample $m$ | 2000 | 2000 | 2000 | 2000 | 2000 | 2000 |  |
In-sample $n$ | 2018 | 2003 | 2081 | 2066 | 2033 | 1986 |  |

_Note_ :Based on average loss and rank, box indicates the favoured model,
dashed box indicates the 2nd ranked model.

##  6 Conclusion

In this paper, an innovative threshold measurement equation is proposed and
its validity is evaluated in a Realized-GARCH framework. Through incorporating
a threshold regression specification, the proposed measurement equation is
capable of capturing the leverage effect in a manner that is different to the
original Realized-GARCH in Hansen et al. (2012). The contemporaneous
dependence between the observed realized measures and hidden volatility is
also successfully modelled in the proposed framework, in a way that is
potentially more flexible compared to the one in Realized-GARCH. This
threshold measurement equation can be employed in various extensions of the
Realized-GARCH, such as Realized-EGARCH, Realized-HAR-GARCH, Realized-
Threshold-GARCH and Realized-CARE, etc.

The estimation of the proposed Realized-T-M-GARCH model employs an adaptive
Bayesian MCMC method, the validity of which is evaluated via a simulation
study. In a empirical study with six market indices and two out-of-sample
sizes, the effectiveness of the proposed model is evaluated. The Realized-T-M-
GARCH model produces $\xi_{1}$ and $\xi_{2}$ estimates that are capable of
adjusting the RV bias dependent on the sign of the lagged return. In addition,
the $\varphi_{1}$ and $\varphi_{2}$ parameters in the threshold measurement
equation present different behaviours which could potentially model the
relationship between the realized measure and volatility in a more flexible
way, compared to the Realized-GARCH. Both observations show that the proposed
threshold measurement equation is capable of capturing the leverage effect
successfully.

The 1% and 2.5% VaR and ES forecasting results lend further evidence on the
usefulness of the proposed framework. Comparing to the Realized-GARCH, the
Realized-T-M-GARCH model produces competitive quantile loss and joint loss
values for an eight year out-of-sample period. Focusing on the high volatility
2008 GFC period, the Realized-T-M-GARCH is in general favoured.

This work could be extended by a number of ways. First, the effectiveness of
incorporating the proposed framework in Realized-EGARCH, Realized-HAR-GARCH,
Realized-Threshold-GARCH and Realized-CARE, etc, could be evaluated and
compared. Second, multiple realized measures could be considered as input to
the proposed threshold measurement equation. The impact of such extended
framework on volatility and tail risk forecasting accuracy could be
investigated. Third, the current threshold specification could be potentially
extended by a smooth transition framework, see the smooth transition GARCH of
González-Rivera (1998) and Anderson et al. (1999), or the smooth transition
dynamic range models of Lin et al. (2012) as examples.

## References

  * Acerbi and Tasche (2002) Acerbi, C. and D. Tasche (2002).  Expected shortfall: A natural coherent alternative to value at risk.  Economic Notes 31(2), 379–388.
  * Andersen and Bollerslev (1998) Andersen, T. G. and T. Bollerslev (1998).  Answering the skeptics: Yes, standard volatility models do provide accurate forecasts.  International economic review, 885–905.
  * Andersen et al. (2003) Andersen, T. G., T. Bollerslev, F. X. Diebold, and P. Labys (2003).  Modeling and forecasting realized volatility.  Econometrica 71(2), 579–625.
  * Anderson et al. (1999) Anderson, H. M., K. Nam, and F. Vahid (1999).  Asymmetric nonlinear smooth transition garch models.  In Nonlinear time series analysis of economic and financial data, pp. 191–207. Springer.
  * Artzner (1997) Artzner, P. (1997).  Thinking coherently.  Risk, 68–71.
  * Artzner et al. (1999) Artzner, P., F. Delbaen, J. Eber, and D. Heath (1999).  Coherent measures of risk.  Mathematical Finance 9(3), 203–228.
  * Black (1976) Black, F. (1976).  Studies of stock market volatility changes.  1976 Proceedings of the American Statistical Association Bisiness and Economic Statistics Section.
  * Bollerslev (1986) Bollerslev, T. (1986).  Generalized autoregressive conditional heteroskedasticity.  Journal of econometrics 31(3), 307–327.
  * Chen et al. (2008) Chen, C. W., R. Gerlach, and E. M. Lin (2008).  Volatility forecasting using threshold heteroskedastic models of the intra-day range.  Computational Statistics & Data Analysis 52(6), 2990–3010.
  * Chen and Watanabe (2019) Chen, C. W. and T. Watanabe (2019).  Bayesian modeling and forecasting of value-at-risk via threshold realized volatility.  Applied Stochastic Models in Business and Industry 35(3), 747–765.
  * Chen et al. (2017) Chen, W. Y., G. W. Peters, R. H. Gerlach, and S. A. Sisson (2017).  Dynamic quantile function models.  arXiv preprint arXiv:1707.02587.
  * Contino and Gerlach (2017) Contino, C. and R. H. Gerlach (2017).  Bayesian tail-risk forecasting using realized garch.  Applied Stochastic Models in Business and Industry 33(2), 213–236.
  * Engle (1982) Engle, R. F. (1982).  Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation.  Econometrica: Journal of the econometric society, 987–1007.
  * Fernández and Steel (1998) Fernández, C. and M. F. Steel (1998).  On bayesian modeling of fat tails and skewness.  Journal of the american statistical association 93(441), 359–371.
  * Fissler and Ziegel (2016) Fissler, T. and J. F. Ziegel (2016).  Higher order elicitability and Osband’s principle.  The Annals of Statistics 44(4), 1680–1707.
  * Gelman et al. (1997) Gelman, A., W. R. Gilks, and G. O. Roberts (1997).  Weak convergence and optimal scaling of random walk metropolis algorithms.  The annals of applied probability 7(1), 110–120.
  * Gerlach and Wang (2016) Gerlach, R. and C. Wang (2016).  Forecasting risk via realized garch, incorporating the realized range.  Quantitative Finance 16(4), 501–511.
  * Gerlach and Wang (2020) Gerlach, R. and C. Wang (2020).  Bayesian semi-parametric realized conditional autoregressive expectile models for tail risk forecasting.  Journal of Financial Econometrics (In Press).
  * Glosten et al. (1993) Glosten, L. R., R. Jagannathan, and D. E. Runkle (1993).  On the relation between the expected value and the volatility of the nominal excess return on stocks.  The journal of finance 48(5), 1779–1801.
  * Gneiting (2011) Gneiting, T. (2011).  Making and evaluating point forecasts.  Journal of the American Statistical Association 106(494), 746–762.
  * Goldfeld and Quandat (1972) Goldfeld, S. M. and R. E. Quandat (1972).  Nonlinear methods in econometrics.
  * González-Rivera (1998) González-Rivera, G. (1998).  Smooth-transition garch models.  Studies in Nonlinear Dynamics & Econometrics 3(2).
  * Granger et al. (1993) Granger, C. W., T. Terasvirta, et al. (1993).  Modelling non-linear economic relationships.  OUP Catalogue.
  * Hansen and Huang (2016) Hansen, P. R. and Z. Huang (2016).  Exponential garch modeling with realized measures of volatility.  Journal of Business & Economic Statistics 34(2), 269–287.
  * Hansen et al. (2012) Hansen, P. R., Z. Huang, and H. H. Shek (2012).  Realized garch: a joint model for returns and realized measures of volatility.  Journal of Applied Econometrics 27(6), 877–906.
  * Huang et al. (2016) Huang, Z., H. Liu, and T. Wang (2016).  Modeling long memory volatility using realized measures of volatility: A realized har garch model.  Economic Modelling 52, 812–821.
  * Li and Li (1996) Li, C. and W. K. Li (1996).  On a double-threshold autoregressive heteroscedastic time series model.  Journal of applied econometrics 11(3), 253–274.
  * Lin et al. (2012) Lin, E. M., C. W. Chen, and R. Gerlach (2012).  Forecasting volatility with asymmetric smooth transition dynamic range models.  International Journal of Forecasting 28(2), 384–399.
  * Nelson (1991) Nelson, D. B. (1991).  Conditional heteroskedasticity in asset returns: A new approach.  Econometrica: Journal of the Econometric Society, 347–370.
  * Poon and Granger (2003) Poon, S.-H. and C. W. Granger (2003).  Forecasting volatility in financial markets: A review.  Journal of economic literature 41(2), 478–539.
  * Silvapulle and Sen (2005) Silvapulle, M. J. and P. K. Sen (2005).  Constrained statistical inference: Inequality, order and shape restrictions.  John Wiley & Sons.
  * Taylor (2019) Taylor, J. W. (2019).  Forecasting value at risk and expected shortfall using a semiparametric approach based on the asymmetric laplace distribution.  Journal of Business & Economic Statistics 37(1), 121–133.
  * Tong (1978) Tong, H. (1978).  On a threshold model.
  * Tong (1990) Tong, H. (1990).  Non-linear time series: a dynamical system approach.  Oxford University Press.
  * Watanabe (2012) Watanabe, T. (2012).  Quantile forecasts of financial returns using realized garch models.  The Japanese Economic Review 63(1), 68–80.


Present address: ]National Institute of Standards and Technology, Boulder,
Colorado 80305 USA Present address: ]Physikalisch Technische Bundesanstalt,
38116 Braunschweig, Germany

# Anomalous broadening in driven dissipative Rydberg systems

E. A. Goldschmidt  United States Army Research Laboratory, Adelphi, Maryland
20783 USA  T. Boulier  Joint Quantum Institute, National Institute of
Standards and Technology and the University of Maryland, Gaithersburg,
Maryland 20899 USA  R. C. Brown  [  Joint Quantum Institute, National
Institute of Standards and Technology and the University of Maryland,
Gaithersburg, Maryland 20899 USA  S. B. Koller  [  Joint Quantum Institute,
National Institute of Standards and Technology and the University of Maryland,
Gaithersburg, Maryland 20899 USA  J. T. Young  Joint Quantum Institute,
National Institute of Standards and Technology and the University of Maryland,
Gaithersburg, Maryland 20899 USA  A. V. Gorshkov  Joint Quantum Institute,
National Institute of Standards and Technology and the University of Maryland,
Gaithersburg, Maryland 20899 USA  Joint Center for Quantum Information and
Computer Science, National Institute of Standards and Technology and the
University of Maryland, College Park, Maryland 20742 USA  S. L. Rolston  Joint
Quantum Institute, National Institute of Standards and Technology and the
University of Maryland, Gaithersburg, Maryland 20899 USA  J. V. Porto
<EMAIL_ADDRESS>Joint Quantum Institute, National Institute of Standards and
Technology and the University of Maryland, Gaithersburg, Maryland 20899 USA

###### Abstract

We observe interaction-induced broadening of the two-photon $5s$-$18s$
transition in 87Rb atoms trapped in a 3D optical lattice. The measured
linewidth increases by nearly two orders of magnitude with increasing atomic
density and excitation strength, with corresponding suppression of resonant
scattering and enhancement of off-resonant scattering. We attribute the
increased linewidth to resonant dipole-dipole interactions of $18s$ atoms with
blackbody induced population in nearby $np$ states. Over a range of initial
atomic densities and excitation strengths, the transition width is described
by a single function of the steady-state density of Rydberg atoms, and the
observed resonant excitation rate corresponds to that of a two-level system
with the measured, rather than natural, linewidth. The broadening mechanism
observed here is likely to have negative implications for many proposals with
coherently interacting Rydberg atoms.

Due to their strong, long-range, coherently-controllable interactions, Rydberg
atoms have been proposed as a basis for quantum information processing and
simulation of many-body physics Weimer  _et al._ (2010); Lee  _et al._ (2012);
Pohl  _et al._ (2011); Saffman  _et al._ (2010). Using the coherent dynamics
of such highly excited atomic states, however, requires addressing challenges
posed by the dense spectrum of Rydberg levels, the detrimental effects of
spontaneous emission, and strong interactions. One approach is to operate on
timescales much faster than the long Rydberg lifetime, typically microseconds
to milliseconds Heidemann  _et al._ (2007); Barredo  _et al._ (2015);
Bendkowsky  _et al._ (2009); Urban  _et al._ (2009). Another proposed approach
is to off-resonantly couple the ground and Rydberg state, admixing a small
amount of the strongly interacting character into the ground state while
substantially reducing spontaneous emission Johnson and Rolston (2010); Henkel
_et al._ (2010); Pupillo  _et al._ (2010); Honer  _et al._ (2010); Glaetzle
_et al._ (2012); van Bijnen and Pohl (2015); Glaetzle  _et al._ (2015);
Bouchoule and Mølmer (2002); Lee  _et al._ (2013); Dauphin  _et al._ (2012);
Glaetzle  _et al._ (2014). This Rydberg-dressed atom approach has been
recently demonstrated with pairs of atoms Jau  _et al._ (2015), but has been
difficult to realize in a many-body context Balewski  _et al._ (2014). A full
understanding of the scope and limitations of these proposals requires
including the effects of spontaneous decay within the dense energy level
structure, which typically cannot be described by a mean-field treatment in
interacting gases due to correlated quantum coherent and dissipative effects.

We study the effect of interactions in a driven, dissipative system of Rydberg
atoms in a 3D optical lattice. We observe significant deviation from the
expected excitation rates both on and off resonance that cannot be explained
by van der Waals interactions or a mean-field treatment of the system. We
attribute these effects to blackbody induced transitions to nearby Rydberg
states of opposite parity, which have large, resonant dipole-dipole
interactions with the state of interest. These off-diagonal exchange
interactions result in complex many-body states of the system. Previous work
has explored the impact of similar, controlled, interactions Anderson  _et
al._ (1998); Park  _et al._ (2011); Günter  _et al._ (2013); Brekke  _et al._
(2012), however the uncontrolled creation of strongly interacting Rydberg
levels due to spontaneous or blackbody processes is typically ignored in
discussions of coherent Rydberg dynamics. These interactions may significantly
modify the parameter regimes available for many-body Rydberg-based systems. In
particular, we show that even at low densities of Rydberg atoms, uncontrolled
production of atoms in other states significantly modifies the energy levels
of the remaining atoms.

Figure 1: a) Level diagram for two-photon excitation to the $18s$ Rydberg
level. Fractional density, $f$, in the $\left|F=2,m_{F}=-2\right\rangle$
ground state is controlled by microwave transfer, remaining atoms are shelved
in the non-participating $\left|F=2,m_{F}=+2\right\rangle$ state (spectrally
resolvable due to a $0.3~{}\rm{mT}$ magnetic field along the optical axis).
Decay from the $18s$ state occurs via many channels, including via Rydberg
$np$ levels and the $5p_{3/2}$ state. Atoms are optically pumped to the
$\left|F=1,2;m_{F}=-1,0\right\rangle$ ground states. b) Example $18s$ spectra
measured as the population in the $m_{F}=0,\pm 1$ states vs two-photon
detuning $\delta$. Blue is $\Omega/2\pi$=3 kHz, $f$=0.3 and red is
$\Omega/2\pi$=140 kHz, $f$=0.75.

We use a state in 87Rb with principal quantum number $n=18$, and relatively
short natural lifetime, $1/\Gamma_{0}=3.5~{}\mu$s (including blackbody
transitions), to study dissipative effects in a Rydberg system. Atoms in the
same Rydberg level interact primarily via the $C_{6}/r^{6}$ van der Waals
interaction, which for the $18s$ state is repulsive and equals the $18s$
linewidth at 800 nm separation. On the other hand, atoms in states of opposite
parity interact via the much larger resonant dipole-dipole interaction,
$\propto C_{3}/r^{3}$, whose angular dependence allows it to be positive or
negative Beterov  _et al._ (2009); Afrousheh  _et al._ (2004) and, for $18s$
interacting with $17p$ or $18p$, equals the $18s$ linewidth at $16~{}\mu$m
separation. Blackbody-induced transitions to other Rydberg levels constitute
$\gtrsim 20~{}\%$ of the decay. We note that molecular resonances can be
ignored due to the low principal quantum number that allows molecule formation
only at extremely high densities Derevianko  _et al._ (2015).

We excite the $18s_{1/2}$ state using a two-photon transition via the 5p1/2
state (Fig. 1a), with intermediate state detuning $\Delta/2\pi\approx 235$ MHz
and independently calibrated single-photon Rabi frequencies,
$\Omega_{1}/2\pi<10$ MHz and $\Omega_{2}/2\pi\approx 7$ MHz sup . The two
excitation lasers are stabilized to the same high-finesse optical cavity with
$<10$ kHz linewidth, and are polarized and tuned to couple the ground
$\left|5s_{1/2};~{}F=2,m_{F}=-2\right\rangle$ hyperfine state to the
$\left|18s_{1/2};~{}F=2,m_{F}=-2\right\rangle$ state with two-photon detuning
$\delta$ and Rabi frequency $\Omega=\Omega_{1}\Omega_{2}/2\Delta$.

The atomic system consists of a Bose-Einstein condensate of $\approx 4\times
10^{4}$ atoms initially in the $\left|F=1,m_{F}=-1\right\rangle$ ground state,
loaded into a 3D optical lattice sup ; Sebby-Strabley  _et al._ (2006). The
lattice provides a minimum separation of 406 nm and additionally suppresses
superradiant Rayleigh scattering on the $5s$-$5p$ transition Inouye  _et al._
(1999). We control the atomic density available for Rydberg excitation with
microwave rapid adiabatic passage that puts a fraction, $f$, of the atoms in
the $\left|F=2,m_{F}=-2\right\rangle$ hyperfine state and shelves all
remaining atoms in the non-participating $\left|F=2,m_{F}=+2\right\rangle$
state. This technique varies the average density, $\rho_{g}=f\times
57~{}\mu$m-3, without altering the geometry of the cloud. We quantify
excitation to the Rydberg state by measuring the population remaining in the
initial state (or, equivalently, pumped into the initially empty $m_{F}=0,\pm
1$ states) following excitation. The ground hyperfine populations are
separated in time-of-flight with a Stern-Gerlach magnetic field gradient and
measured via absorption imaging.

Figure 2: a) Measured optical pumping rate for $5s$-$5p$ field only (red) vs
$\Omega_{1}$ (top axis) and two-photon resonant excitation with lower rate
subtracted (black) vs $\Omega$ (bottom axis). Black line is a linear fit to
$R_{0}$ and red line is the calculated $5s$-$5p$ rate with no adjustable
parameters. Gray dashed line is the expected single particle rate. b) Measured
width $\Gamma$ in units of the natural linewidth $\Gamma_{0}=2\pi\times
45~{}\rm{kHz}$ vs $\Omega$ for different fractional densities. Dashed line is
linear scaling. Error bars represent statistical fitting uncertainties.

We measure the Rydberg excitation rate $R$ (proportional to the optically
pumped fraction following a fixed excitation time at least several $18s$
lifetimes, but shorter than the time to depump the initial state) as a
function of $\delta$, $\Omega$ and $\rho_{g}$. Observed lineshapes are
symmetric and well-characterized by Lorentzians (Fig. 1b)

$R=\frac{R_{0}}{1+4\delta^{2}/\Gamma^{2}}.$ (1)

We fit a decaying exponential to the population remaining in the initial state
as a function of excitation time for two-photon excitation with $\delta=0$ and
for the lower $5s$-$5p$ field alone (Fig. 2a). For each $\Omega$ and
$\rho_{g}$, we extract the resonant excitation rate $R_{0}$ by subtracting the
$5s$-$5p$ optical pumping rate from the measured two-photon rate and scaling
by the $45~{}\%$ fraction that decays to states other than the initial state.
The linewidth $\Gamma$ is determined from a Lorentzian fit to the optical
pumping as a function of $\delta$ (Fig. 2b). We observe that $\Gamma$
increases dramatically with both $\Omega$ and $\rho_{g}$, reaching values as
large as $\approx 200\Gamma_{0}$. At small $\Omega$ and $\rho_{g}$ the
narrowest observed linewidth is $\approx 3\Gamma_{0}$ and the residual
broadening is attributed to inhomogeneities such as optical trapping light
shifts and laser frequency noise.

Remarkably, $R_{0}$ is linear in $\Omega$ (with slope that depends on
$\rho_{g}$) and shows no sign of saturation up to $\Omega=3\Gamma_{0}$ (Fig.
2a). This behavior is inconsistent with standard single-atom theory and purely
inhomogeneous broadening, which predicts faster excitation that depends on
$\Omega^{2}$ for small $\Omega$ and saturates at large $\Omega$ (Fig. 2a
dashed line). The observed $R_{0}$ corresponds to a single-atom theory
assuming the measured $\Gamma$ as the transition linewidth:
$R_{0}\approx\Omega^{2}/\Gamma$.

To determine whether the observed broadening corresponds to a concomitant
shortening of the $18s$ lifetime, as would broadening due to superradiance
Wang  _et al._ (2007), we collect fluorescence emitted on the
$5p_{3/2}$-$5s_{1/2}$ transition (Fig. 1a). The fluorescence, which scales
with the optical pumping signal and is proportional to the number of $18s$
atoms, is collected by a lens relay system (NA=0.12) with an interference
filter to block the $5s$-$5p_{1/2}$ excitation light, detected by a single
photon avalanche diode, time-tagged with 21-ns resolution, and summed over
many excitation pulses. The observed lifetime, measured as the decay in
detected photons after extinguishing the excitation light, is consistent with
the 3.5 $\mu$s natural lifetime and independent of $\Omega$ (see Fig. 4 inset
and sup  for more information). This result is consistent with previous
observations of the suppression of superradiance due to driven dipole
interactions Pellegrino  _et al._ (2014). The confirmation of the natural
lifetime, along with the lack of saturation of the optical pumping, rules out
superradiance and suggests the broadening is due to rapid dephasing of the
optical coherence. In addition, confirmation of the lifetime allows an
estimate of the steady-state $18s$ population.

The steady-state density of $18s$ atoms, under resonant excitation, is the
atomic density, $\rho_{g}$, scaled by the ratio of the excitation rate,
$R_{0}$, to the decay rate, $\Gamma_{0}$:
$\rho_{18s}=\rho_{g}R_{0}/\Gamma_{0}$. The steady-state densities in nearby
$np$ states are equal to $\rho_{18s}$ scaled by the ratios of the $18s-np$
transition rates to the $np$ decay rates:
$\rho_{np}=\rho_{18s}b_{np}\Gamma_{0}/\Gamma_{np}=\rho_{g}R_{0}b_{np}/\Gamma_{np}$,
where $b_{np}$ are the branching ratios from $18s$ to $np$ (dominated by
blackbody transitions to $17p$ and $18p$ states), and $\Gamma_{np}$ are the
decay rates of the $np$ states (including blackbody transitions).

We observe $\Gamma$ as large as 8 MHz, inconsistent with the 1.9 MHz van der
Waals shift expected at our highest $18s$ densities sup . In addition, the
observed lineshapes are symmetric, inconsistent with the repulsive van der
Waals interaction. Also, the broadening depends only on the average density
$\rho_{g}$ and is independent of the microscopic configuration, which we alter
by transferring atoms in every other lattice site in 2D to the shelving state
Lee  _et al._ (2007). The width, lineshape symmetry, and insensitivity to
nearest-neighbor spacing are consistent with the larger, longer-range,
symmetric, dipole-dipole interaction between states of opposite parity.

Figure 3: (a) Measured width $\Gamma$ vs $\Omega\sqrt{\rho_{g}\beta_{3}}$ and
(b) resonant excitation rate $R_{0}$ vs $\Omega/\sqrt{\rho_{g}\beta_{3}}$ at
different two-photon Rabi frequencies and densities. Linear scalings with unit
slope are indicated with solid lines. Error bars represent statistical fitting
uncertainties.

For broadening due to dipole interactions, we expect a width of order
$\sum\left|C_{3}^{(np)}\right|\rho_{np}$, where the sum is over the $np$
states, which have different effective interaction strengths and branching
ratios. This can be rewritten as $\beta_{3}\rho_{g}R_{0}$ using the expression
above for $\rho_{np}$ and defining the quantity
$\beta_{3}=\sum\left|C_{3}^{(np)}\right|b_{np}/\Gamma_{np}=116~{}\mu\mathrm{m}^{3}$
(including the root-mean-squared average of the angular dependence of
$C_{3}$). Combined with the observed relation $R_{0}\approx\Omega^{2}/\Gamma$,
this provides expressions for $\Gamma$ and $R_{0}$ in terms of independently
controlled variables $\Omega$ and $\rho_{g}$:

$\begin{split}\Gamma\approx\Omega\sqrt{\rho_{g}\beta_{3}},\\\
R_{0}\approx\frac{\Omega}{\sqrt{\rho_{g}\beta_{3}}}.\end{split}$ (2)

$\Gamma$ and $R_{0}$ are plotted in terms of these expressions in Fig. 3. The
data not only collapse to approximately linear curves over two orders of
magnitude, but the magnitude is well described by the dipole-dipole energy
scale characterized by the independently calculated factor $\beta_{3}$.
(Neither of these features is present for scaling with the van der Waals
interaction sup ). This agreement is highly suggestive of a broadening
mechanism dominated by dipole interactions with contaminant states. The
fluctuating microscopic configuration of $np$ states being populated and
decaying leads to dephasing that is not accompanied by either saturated
optical pumping or shortening of the lifetime.

This broadening mechanism requires some initial time to populate the
contaminant states. We study the time dynamics of resonant and detuned
excitation using the fluorescence on the $5p_{3/2}$-$5s_{1/2}$ transition.
Fig. 4 shows the fluorescence, converted into a number of $18s$ atoms, as a
function of time for excitation at different detunings with $\Omega/2\pi=140$
kHz and $f=1$. At the two non-zero detunings, the population reaches a
significant fraction of the resonantly excited population in a few $\mu$s.
This is in stark contrast to the expected single-atom scattering times of
$3~{}\mathrm{ms}$ and $11~{}\mathrm{ms}$ for these nominally far detuned cases
($\tau_{s}=4\delta^{2}/\Gamma_{0}\Omega^{2}$ for
$\delta\gg\Omega,\Gamma_{0}$). The faster excitation off resonance leads to
observed $18s$ populations larger by a factor of $\gtrsim 30$ than expected.
For the resonant case, on the other hand, the observed population is smaller
than expected from a single-atom picture by a factor of $\gtrsim 10$, which
cannot be explained by van der Waals interactions alone. The expected non-
interacting excitation rates, both on and off resonance, are central to the
feasibility of both dressed Rydberg proposals and Rydberg quantum gate
implementations. We observe substantial deviations from these expected rates,
which must be addressed for a full analysis of any Rydberg system.

Figure 4: Average fluorescence photon counts (converted to number of $18s$
atoms) as the excitation turns on for different two-photon detunings $\delta$
for $\Omega/2\pi$=140 kHz and $f=1$. Dashed lines indicate the number of
excited atoms expected for the non-interacting, single particle scattering
rates at those detunings. Arrows indicate the expected times to create the
first contaminant atom for the detuned cases (see text). Error bars represent
shot noise. Inset: Decay following short ($2~{}\mu$s) pulses for different
two-photon Rabi frequencies with natural lifetime plotted for reference. All
regimes are consistent with the $18s$ natural lifetime.

This shortening of the off-resonant scattering time may not be a problem in
few atom systems such as arrays of microtraps or other 2D systems Jau  _et
al._ (2015); Zeiher  _et al._ (2015); Xia  _et al._ (2015), but it is likely
problematic for implementing Rydberg-dressed atom proposals in large, many-
body systems Johnson and Rolston (2010); Henkel  _et al._ (2010); Pupillo  _et
al._ (2010); Honer  _et al._ (2010). The time until the creation of the first
contaminant atom is $\tau=\tau_{s}/bN_{0}$, where $b\gtrsim 20~{}\%$ is the
branching to contaminant states and $N_{0}$ is the total number of atoms.
Interaction with the first $np$ atom allows excitation of other atoms at a
faster, resonant rate, leading to additional $np$ atoms that, in turn,
increase the number of atoms resonantly excited, similar to Rydberg
aggregation at shorter timescales due to van der Waals shifts Schempp  _et
al._ (2014); Urvoy  _et al._ (2015); Malossi  _et al._ (2014). The long-range
nature of the dipole interaction causes aggregation on length scales
comparable to typical experimental system sizes, leading to rapid broadening
over the entire ensemble. This simple timescale estimate gives a qualitative
understanding of the early time dynamics, and future work to develop a full
microscopic model will hopefully provide quantitative descriptions of both the
dynamic and steady-state behavior sup . Finally, the magnitude of the
uncontrolled interactions with contaminant atoms is large compared to the
interactions in a Rydberg-dressed approach. In particular, a dressed atom’s
uncontrolled dipole interaction with a contaminant atom is larger than its
interaction with another dressed atom for $\rho_{np}>\Omega^{2}/\delta C_{3}$
Balewski  _et al._ (2014), which is quickly exceeded under reasonable
experimental conditions sup .

In conclusion, we report experimental observation of large spectral broadening
of a Rydberg transition modifying the scattering rate both on and off
resonance. We infer this effect results from the uncontrolled buildup of atoms
in nearby Rydberg states. Resonant dipole-dipole interaction with those states
causes dephasing and broadens the driven transition. Any single-atom approach
to this problem is inherently nonlinear, as the broadening depends on the
excited population, leading to distinctly non-Lorentzian lineshapes that
contradict our observations. Mean-field approaches fail because the off-
diagonal interaction requires single-atom coherences between the driven state
and contaminant states, which do not develop under blackbody induced
population of the contaminant states sup . This suggests the importance of
correlations and is the focus of future theoretical efforts. Nonetheless,
independent of a microscopic model, a simple analysis supported by
experimental observation suggests the time available for coherent
manipulations is much shorter than the expected single-atom scattering time,
placing significant constraints on Rydberg dressing proposals. Importantly,
the mechanism described here scales unfavorably with principal quantum number
sup  and implies the need to account for even a small number of impurity
Rydberg atoms when considering interactions in dense gases. And although we
have focused on exciting to an s-state with contaminant p-states, this
mechanism is similarly present for excitation to any Rydberg state, which will
populate nearby states of opposite parity. We note similar broadening has been
observed in Rydberg transitions in strontium DeSalvo  _et al._ (2015).

###### Acknowledgements.

The authors thank T. C. Killian, F. B. Dunning, A. Browaeys, M. Foss-Feig, A.
Hu, R. M. Wilson, and Z. X. Gong for helpful discussions. This work was
partially supported by NSF PIF, AFOSR, ARO, ARL-CDQI, and NSF PFC at JQI.

## References

  * Weimer  _et al._ (2010) H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, Nature Physics 6, 382 (2010).
  * Lee  _et al._ (2012) T. E. Lee, H. Häffner, and M. C. Cross, Physical Review Letters 108, 023602 (2012).
  * Pohl  _et al._ (2011) T. Pohl, C. S. Adams, and H. R. Sadephpour, Journal of Physics B: Atomic, Molecular and Optical Physics 44, 180201 (2011).
  * Saffman  _et al._ (2010) M. Saffman, T. G. Walker, and K. M$\o$lmer, Reviews of Modern Physics 82, 2313 (2010).
  * Heidemann  _et al._ (2007) R. Heidemann, U. Raitzsch, V. Bendkowsky, B. Butscher, R. Löw, L. Santos, and T. Pfau, Physical Review Letters 99, 163601 (2007).
  * Barredo  _et al._ (2015) D. Barredo, H. Labuhn, S. Ravets, T. Lahaye, A. Browaeys, and C. S. Adams, Physical Review Letters 114, 113002 (2015).
  * Bendkowsky  _et al._ (2009) V. Bendkowsky, B. Butscher, J. Nipper, J. P. Shaffer, R. Löw, and T. Pfau, Nature 458, 1005 (2009).
  * Urban  _et al._ (2009) E. Urban, T. A. Johnson, T. Henage, L. Isenhower, D. Yavuz, T. Walker, and M. Saffman, Nature Physics 5, 110 (2009).
  * Johnson and Rolston (2010) J. E. Johnson and S. L. Rolston, Physical Review A 82, 033412 (2010).
  * Henkel  _et al._ (2010) N. Henkel, R. Nath, and T. Pohl, Physical Review Letters 104, 195302 (2010).
  * Pupillo  _et al._ (2010) G. Pupillo, A. Micheli, M. Boninsegni, I. Lesanovsky, and P. Zoller, Physical Review Letters 104, 223002 (2010).
  * Honer  _et al._ (2010) J. Honer, H. Weimer, T. Pfau, and H. P. Büchler, Physical Review Letters 105, 160404 (2010).
  * Glaetzle  _et al._ (2012) A. W. Glaetzle, R. Nath, B. Zhao, G. Pupillo, and P. Zoller, Physical Review A 86, 043403 (2012).
  * van Bijnen and Pohl (2015) R. M. W. van Bijnen and T. Pohl, Physical Review Letters 114, 243002 (2015).
  * Glaetzle  _et al._ (2015) A. W. Glaetzle, M. Dalmonte, R. Nath, C. Gross, I. Bloch, and P. Zoller, Physical Review Letters 114, 173002 (2015).
  * Bouchoule and Mølmer (2002) I. Bouchoule and K. Mølmer, Physical Review A 65, 041803 (2002).
  * Lee  _et al._ (2013) T. E. Lee, S. Gopalakrishnan, and M. D. Lukin, Physical Review Letters 110, 257204 (2013).
  * Dauphin  _et al._ (2012) A. Dauphin, M. Müller, and M. A. Martin-Delgado, Physical Review A 86, 053618 (2012).
  * Glaetzle  _et al._ (2014) A. W. Glaetzle, M. Dalmonte, R. Nath, I. Rousochatzakis, R. Moessner, and P. Zoller, Physical Review X 4, 041037 (2014).
  * Jau  _et al._ (2015) Y.-Y. Jau, A. M. Hankin, T. Keating, I. H. Deutsch, and G. W. Biedermann, Nature Physics 12, 71 (2016) .
  * Balewski  _et al._ (2014) J. B. Balewski, A. T. Krupp, A. Gaj, S. Hofferberth, R. Löw, and T. Pfau, New Journal of Physics 16, 063012 (2014).
  * Anderson  _et al._ (1998) W. R. Anderson, J. R. Veale, and T. F. Gallagher, Physical Review Letters 80, 249 (1998).
  * Park  _et al._ (2011) H. Park, P. J. Tanner, B. J. Claessens, E. S. Shuman, and T. F. Gallagher, Physical Review A 84, 022704 (2011).
  * Günter  _et al._ (2013) G. Günter, H. Schempp, M. Robert-de Saint-Vincent, V. Gavryusev, S. Helmrich, C. S. Hofmann, S. Whitlock, and M. Weidemüller, Science 342, 954 (2013).
  * Brekke  _et al._ (2012) E. Brekke, J. O. Day, and T. G. Walker, Phys. Rev. A 86, 033406 (2012).
  * Beterov  _et al._ (2009) I. I. Beterov, I. I. Ryabtsev, D. B. Tretyakov, and V. M. Entin, Physical Review A 79, 052504 (2009).
  * Afrousheh  _et al._ (2004) K. Afrousheh, P. Bohlouli-Zanjani, D. Vagale, A. Mugford, M. Fedorov, and J. D. D. Martin, Physical Review Letters 93, 233001 (2004).
  * Derevianko  _et al._ (2015) A. Derevianko, P. Kómár, T. Topcu, R. M. Kroeze, and M. D. Lukin, Physical Review A 92, 063419 (2015) .
  * (29) See supplementary material.
  * Li  _et al._ (2003) W. Li, I. Mourachko, M. W. Noel, and T. F. Gallagher, Physical Review A 67, 052502 (2003).
  * Sebby-Strabley  _et al._ (2006) J. Sebby-Strabley, M. Anderlini, P. S. Jessen, and J. V. Porto, Physical Review A 73, 033605 (2006).
  * Inouye  _et al._ (1999) S. Inouye, A. Chikkatur, D. Stamper-Kurn, J. Stenger, D. Pritchard, and W. Ketterle, Science 285, 571 (1999).
  * Wang  _et al._ (2007) T. Wang, S. F. Yelin, R. Côté, E. E. Eyler, S. M. Farooqi, P. L. Gould, M. Koštrun, D. Tong, and D. Vrinceanu, Phys. Rev. A 75, 033802 (2007).
  * Pellegrino  _et al._ (2014) J. Pellegrino, R. Bourgain, S. Jennewein, Y. R. P. Sortais, A. Browaeys, S. D. Jenkins, and J. Ruostekoski, Phys. Rev. Lett. 113, 133602 (2014).
  * Lee  _et al._ (2007) P. J. Lee, M. Anderlini, B. L. Brown, J. Sebby-Strabley, W. D. Phillips, and J. V. Porto, Physical Review Letters 99, 020402 (2007).
  * Zeiher  _et al._ (2015) J. Zeiher, P. Schauß, S. Hild, T. Macrì, I. Bloch, and C. Gross, Physical Review X 5, 031015 (2015).
  * Xia  _et al._ (2015) T. Xia, M. Lichtman, K. Maller, A. W. Carr, M. J. Piotrowicz, L. Isenhower, and M. Saffman, Physical Review Letters 114, 100503 (2015).
  * Schempp  _et al._ (2014) H. Schempp, G. Günter, M. Robert-de Saint-Vincent, C. S. Hofmann, D. Breyel, A. Komnik, D. W. Schönleber, M. Gärttner, J. Evers, S. Whitlock, and M. Weidemüller, Physical Review Letters 112, 013002 (2014).
  * Urvoy  _et al._ (2015) A. Urvoy, F. Ripka, I. Lesanovsky, D. Booth, J. P. Shaffer, T. Pfau, and R. Löw, Physical Review Letters 114, 203002 (2015).
  * Malossi  _et al._ (2014) N. Malossi, M. M. Valado, S. Scotto, P. Huillery, P. Pillet, D. Ciampini, E. Arimondo, and O. Morsch, Physical Review Letters 113, 023006 (2014).
  * DeSalvo  _et al._ (2015) B. J. DeSalvo, J. A. Aman, C. Gaul, T. Pohl, S. Yoshida, J. Burgdörfer, K. R. A. Hazzard, F. B. Dunning, and T. C. Killian, Physical Review A 93, 022709 (2016).

Supplementary material: Anomalous broadening in driven dissipative Rydberg
systems

##

## Rabi frequency calibration

We calibrate the Rabi frequencies of the two optical fields using the light
shifts due to each field. For the 795-nm laser, which is detuned
$\Delta/2\pi\approx 235$ MHz below the
$\left|5s;F=2,m_{F}=-2\right\rangle\rightarrow\left|5p_{1/2};F=1,m_{F}=-1\right\rangle$
transition, we observe the intensity dependent shift, $\delta_{795}$, of the
$\left|5s;F=2,m_{F}=-2\right\rangle\rightarrow\left|5s;F=1,m_{F}=-1\right\rangle$
microwave transition. The light shift $\delta_{795}=\Omega_{1}^{2}/4\Delta$ is
a factor of $\Delta_{\mathrm{hyperfine}}/\Delta\approx 30$ larger for the
$\left|5s;F=2,m_{F}=-2\right\rangle$ state compared to the
$\left|5s;F=1,m_{F}=-1\right\rangle$ state so we take the effect on the
microwave transition as the light shift on the
$\left|5s;F=2,m_{F}=-2\right\rangle$ state. We observe a shift that is linear
in intensity $I$, as expected, and use this to calibrate the Rabi frequency
$\Omega_{1}(I)=2\sqrt{\delta_{795}(I)\Delta}$ which we vary from $2\pi\times
0.2$ MHz to $2\pi\times 10$ MHz to control the two-photon Rabi frequency.

To obtain the Rabi frequency of the 485-nm field, detuned $\Delta/2\pi\approx
235$ MHz above the
$\left|5p_{1/2};F=1,m_{F}=-1\right\rangle\rightarrow\left|18s;F=2,m_{F}=-2\right\rangle$
transition, we observe the intensity dependent shift, $\delta_{485}$, of the
two-photon
$\left|5s;F=2,m_{F}=-2\right\rangle\rightarrow\left|18s;F=2,m_{F}=-2\right\rangle$
transition keeping the 795-nm light intensity (and frequency) constant. The
485-nm field is far from any transition coupling the ground $5s$ state to any
optically excited state, so the shift of the
$\left|18s;F=2,m_{F}=-2\right\rangle$ state is the dominant contribution and
$\delta_{485}=\Omega_{2}^{2}/4\Delta$. We observe a shift that is linear in
intensity, as expected, and use this to calibrate the Rabi frequency
$\Omega_{2}(I)=2\sqrt{\delta_{485}(I)\Delta}$ which we hold at the maximum
value given available laser power and beam diameter, $\Omega_{2}/2\pi\approx
7$ MHz.

## Lattice filling

Loading $4\times 10^{4}$ atoms into the 3D optical lattice as described in the
main text leads to overfilling (more than one atom per site) in a fraction of
the lattice sites. However, we determined that the filling fraction does not
affect the broadening, which is only determined by the overall density. We
made this determination by comparing a Mott insulator state with no more than
one atom per site to a state with a Poissonian distribution of atoms per site,
but the same total number of atoms in the participating ground state (and thus
the same global density). Only the overall atom number was pertinent to the
broadening. We also compared random transfer of half of the atoms to the
participating ground state with transfer of all atoms on every other site in
2D in a checkerboard fashion Lee  _et al._ (2007). These two cases lead to a
factor of two difference in the filling fraction per lattice site. We again
saw that the total atom number controlled the broadening without regard to the
microscopic configuration.

##  $18s$ Lifetime

We observe 780-nm fluorescence on the $5p_{3/2}-5s$ transition following
excitation to the 18s Rydberg state and fit a decaying exponential to extract
a lifetime. The bulk of the fluorescence is due to the $\approx 30~{}\%$ of
the decay via $18s-5p_{3/2}-5s$ for which the 28-ns lifetime of the $5p_{3/2}$
is negligible compared to $\tau_{0}=3.5~{}\mu\mathrm{s}$, the natural lifetime
of the 18s. However, $\approx 1~{}\%$ of the population decays via channels of
the form $18s-np-n^{\prime}s-5p_{3/2}-5s$ for which the lifetimes of the
intermediate states delay the final 780-nm photon leading to effective
lifetimes $\approx 4\times$ longer Beterov  _et al._ (2009) and an increase in
the measured lifetime by $\approx 10~{}\%$. In addition, radiation trapping,
in which photons are reabsorbed and reemitted one or more times before leaving
the cloud, can alter the measured lifetime. By measuring fluorescence
following resonant excitation on the $5s_{1/2}-5p_{3/2}$ transition, we place
a limit of $<0.5~{}\mu\mathrm{s}$, shorter than the measured lifetimes, on the
characteristic time for radiation trapping, and thus expect little or no
alteration of the measured lifetime due to this effect.

In all cases studied, including resonant excitation at different two-photon
Rabi frequencies and detuned excitation at different detunings, the extracted
lifetime is $8~{}\%$ to $20~{}\%$ longer than $\tau_{0}$, consistent with the
natural lifetime and suggesting no shortening of the lifetime due to
superradiance or other purely homogeneous effects. In addition, the amount of
fluorescence is consistent with the optical pumping signal under all different
conditions. This further suggests that superradiance is not a significant
effect as it would lead to a Rydberg density dependence of the fractional
decay via the $5p_{3/2}$ state compared to other intermediate states.

## Dipole vs van der Waals scaling

In order to calculate $\beta_{3}$, we have only included the $17p$ and $18p$
states and have excluded interactions that do not conserve total magnetic
quantum number. The former should add $\approx 1~{}\%$ to the total
interaction and the latter are not resonant due to Zeeman splitting caused by
a non-zero magnetic field. However some $m$ non-conserving interactions have
an energy mismatch that is less than the largest observed linewidth and may
play a role.

If the broadening were due to van der Waals interactions between $18s$ atoms,
one expects scaling $\Gamma=C_{6}\rho_{18s}^{2}$. We define an interaction
volume $\beta_{6}=\sqrt{C_{6}/\Gamma_{0}}$, in which case we assume
$\Gamma=\beta_{6}^{2}\rho_{g}^{2}R_{0}^{2}/\Gamma_{0}$. Combined with the
observed relation $R_{0}=\Omega^{2}/\Gamma$, we express $\Gamma$ and $R_{0}$
in terms of the independently controlled variables $\rho_{g}$ and $\Omega$:

$\begin{split}&\Gamma=\left(\Omega^{4}/\Gamma_{0}\right)^{1/3}\
\left(\rho_{g}\beta_{6}\right)^{2/3}\\\
&R_{0}=\left(\Omega^{2}\Gamma_{0}\right)^{1/3}\
\left(\rho_{g}\beta_{6}\right)^{-2/3}.\end{split}$ (S1)

Figure S1 shows $\Gamma$ and $R_{0}$ in terms of these expressions. There is a
large mismatch in values between the data and the van der Waals scaling. In
addition the data is not linear in the expressions and the rate $R_{0}$ does
not collapse to a single function for all $\Omega$ and $\rho_{g}$.

Figure S1: (a) Measured $\Gamma$ and (b) $R_{0}$ for different $\Omega$ and
$\rho_{g}$ vs expected scaling for van der Waals dominated interaction. Linear
scalings with unit slope are indicated with solid lines.

## Dipole broadening in various regimes

We consider how broadening due to dipole interaction with spontaneously
created atoms in nearby Rydberg states affects a variety of possible operating
regimes. For systems with small atom number, the time until the first atom is
created in a nearby Rydberg state can be long compared to the lifetime if
$bN_{0}\ll 4\delta^{2}/\Omega^{2}$. However, since the dressed Rydberg
interaction is of order $\Omega^{4}/\delta^{3}$ it may be difficult to access
interesting many-body effects in this regime.

Additionally, this effect remains at higher principal quantum number:
$\beta_{3}\propto n^{*7}$ where $n^{*}=n-\delta_{qd}$ is the effective
principal quantum number and $\delta_{qd}$ is the quantum defect Li  _et al._
(2003). The contaminant density at which the dipole interaction overtakes the
dressed interaction ($\Omega^{2}/\delta C_{3}$) scales as $n^{*-4}$. Assuming
reasonable limits on attainable Rabi frequency ($\Omega/2\pi<100~{}\rm{MHz}$)
and dressed interaction strength ($\Omega^{4}/\delta^{3}<2\pi\times
10~{}\rm{kHz}$), this corresponds to a single spontaneous Rydberg atom in a
typical 10-$\mu\rm{m}$ diameter ultracold atom system for $n>40$. Systems at
higher principal quantum number may allow some additional time in an absolute
sense due to their longer lifetimes (which scale as $<n^{*3}$), but are still
subject to this effect and cannot access the expected single-atom scattering
rates.

Since the van der Waals interaction is derived from the dipole interaction
with nearby states, any configuration with large van der Waals interaction
must have dipole interactions at least as large. Even in the case of operation
at a Förster resonance to maximize the interaction Pupillo  _et al._ (2010);
van Bijnen and Pohl (2015), the uncontrolled, spontaneously created population
in nearby Rydberg states causes interactions larger than any dressed
interaction and must be considered in determining the timescales available for
coherent operations.

## Inhomogeneous Gutzwiller Mean-Field Theory

Treating each atom as a three level system with states $g$, $p$, and $s$,
corresponding to ground, $18s$, and $np$ states respectively, we model the
system with the following Hamiltonian:

$H=-\delta\sum_{i}\sigma_{i}^{ss}+\Omega\sum_{i}(\sigma_{i}^{sg}+\sigma_{i}^{gs})+\sum_{ij}V_{ij}\sigma_{i}^{sp}\sigma_{j}^{ps}+h.c.,$
(S2)

where $V_{ij}=\frac{C_{3}}{r_{ij}^{3}}(1-3\cos^{2}\theta)$ is the dipole-
dipole interaction with $\theta$ the angle between the $z$-axis and the
relative position $\mathbf{r}_{ij}$. The system evolves according to the
master equation

$\dot{\rho}=-i[H,\rho]+\mathcal{L}_{s}+\mathcal{L}_{p}+\mathcal{L}_{R},$ (S3a)
$\mathcal{L}_{s}=\Gamma_{s}\sum_{i}\left(\sigma_{i}^{gs}\rho\sigma_{i}^{sg}-\frac{1}{2}\\{\sigma_{i}^{ss},\rho\\}\right),$
(S3b)
$\mathcal{L}_{p}=\Gamma_{p}\sum_{i}\left(\sigma_{i}^{gp}\rho\sigma_{i}^{pg}-\frac{1}{2}\\{\sigma_{i}^{pp},\rho\\}\right),$
(S3c)
$\mathcal{L}_{R}=\Gamma_{R}\sum_{i}\left(\sigma_{i}^{ps}\rho\sigma_{i}^{sp}-\frac{1}{2}\\{\sigma_{i}^{ss},\rho\\}\right),$
(S3d)

where $\mathcal{L}_{s}$, $\mathcal{L}_{p}$, and $\mathcal{L}_{R}$ are Lindblad
terms corresponding to decay from $s$ to $g$, from $p$ to $g$, and from $s$ to
$p$ respectively.

Using an inhomogeneous Gutzwiller mean-field approximation, we assume the
density matrix has the form

$\rho=\bigotimes_{i}\rho_{i},$ (S4)

which assumes there are no correlations between different atoms. The method is
inhomogeneous in the sense that each atom has its own density matrix, whereas
in homogeneous Gutzwiller mean-field theory all atoms have the same density
matrix. This results in an effective local Hamiltonian

$H_{i}^{eff}=-\delta\sigma_{i}^{ss}+\frac{\Omega}{2}(\sigma_{i}^{sg}+\sigma_{i}^{gs})+\sum_{j}V_{ij}\sigma_{i}^{sp}\langle\sigma_{j}^{ps}\rangle+h.c.$
(S5)

In this picture, the interactions behave as an effective driving term between
the $s$ and $p$ states whose strength and phase are determined by the
$\langle\sigma^{ps}\rangle$ coherences of the surrounding atoms.

We determine the steady state of these equations numerically by initializing a
cubic lattice of randomized density matrices for each site and evolving the
system according to the master equation and effective Hamiltonian. This was
done for a variety of $\Omega$ and $\delta$ based on the experimental
measurements. The interaction strengths were reduced due to computational
constraints, but the nearest neighbor interaction strength remained at least
two orders of magnitude above $\Omega$ and all decay rates, compared to four
in the experiment.

In all cases, we find that the $\langle\sigma^{ps}\rangle$ coherences all
decay to zero in steady state, in which case the system behaves as if there
are no interactions. This is a consequence of the flip-flop interactions and
would not occur if the interactions were of the form

$\sum_{ij}V_{ij}\sigma^{ss}_{i}\sigma^{pp}_{j}.$ (S6)

Collective decay between the $s$ and $p$ states does result in nonzero
$\langle\sigma^{ps}\rangle$ coherences in steady state, but the effect of
interactions in this case is small and the experimental results indicate that
collective decay is not the source of the observed broadening.

Future work to solve this model using techniques beyond Gutzwiller mean-field
will hopefully provide a more full understanding of the observed broadening
and inform a more detailed study of the early time dynamics of the system.