A FORMAL METHOD FOR IDENTIFYING DISTINCT STATES OF VARIABILITY IN TIME-VARYING SOURCES: SGR A* AS AN EXAMPLE

Let us consider a hypothetical flux time series where y_t is the flux measured at time t, and the data are well described by a stochastic model of the form

$$\begin{equation} y_t = c_1 + \phi y_{t-1} + \epsilon _t, \end{equation} \tag{ 1 }$$

where c₁ and |ϕ| < 1 are constants and $\epsilon \sim \mathcal {N} (0,\sigma ^2)$. Let us further assume that we would like to test whether a model that allows a change in mean flux gives a better description of the data. For a permanent change of the mean flux, we could just write down a different parameter c₂ for some t > t₁,

$$\begin{equation} y_t = c_2 + \phi y_{t-1} + \epsilon _t. \end{equation} \tag{ 2 }$$

We are, however, interested in a situation where we can jump back and forth between both models. We will refer to the different models as different states, and at a given time, the source emitting the flux is in state s_t = j, where j = 1 or 2 in our example. We will assume that this state is not directly observable (e.g., with spectroscopic observations) but has to be inferred from the light curve itself.

The unobserved state is modeled using a Markov switching approach. Concretely, let p_ij be a matrix that reflects the transition probabilities of switching to another state or remaining in the current state,

$$\begin{equation} p_{ij} = {\rm Prob}(s_t = j | s_{t-1} = i), \end{equation} \tag{ 3 }$$

and let Π_jt be the probability of being in state j at time t,

$$\begin{equation} \Pi _{jt} = {\rm Prob}(s_t = j | y_t). \end{equation} \tag{ 4 }$$

The probability density of y_t for our model is then

$$\begin{eqnarray} f(y_t | s_t=j,y_{t-1}) &=& \frac{1}{\sqrt{2\pi \sigma ^2}} \exp \left({\frac{-(y_t-c_j-\phi y_{t-1})^2}{2 \sigma ^2}}\right) \nonumber\\ &=& f_{jt}. \end{eqnarray} \tag{ 5 }$$

With these notations, the conditional probability density of y_t can be written as

$$\begin{equation} h(y_t | y_{t-1};\theta) = \sum _{i,j} p_{ij} \Pi _{i t-1} f_{jt}, \end{equation} \tag{ 6 }$$

where θ means the set of parameters in the model (here, c_j, ϕ, and p_ij). Since Π_jt is not directly observable, it has to be recursively calculated by observing that

$$\begin{equation} \Pi _{jt} = \frac{\sum _i p_{ij} \Pi _{it-1} f_{jt}}{h(y_t | y_{t-1}; \theta)}. \end{equation} \tag{ 7 }$$

Starting at t₀, we can execute the iteration to solve for all Π_jt. This means that the final log-likelihood function of the whole light curve is given by

$$\begin{equation} \log h(y_1,y_2,y_3,...,y_T | y_0; \theta) = \sum _{t=1}^T \log h(y_t | y_{t-1};\theta). \end{equation} \tag{ 8 }$$

The maximum of this log-likelihood function gives the preferred values for c_j, ϕ, and p_ij.

While the simple model elucidates the construction of the likelihood function, it is not well suited to be applied to light curves from Sgr A*. The main reason is that the real data set shows unevenly sampled measurements with big gaps representing limited telescope time, the night/day cycle, the observability of the Galactic center from the ground throughout a year, measurements of sky background, and instrument failures. This requires us to use a more appropriate model that is time continuous and not discrete as in the above example. Furthermore, the flux density distribution is not Gaussian, and the data contain measurement noise. This is where the astronomical application of the regime switching approach departs from mathematical finance. We will deal with these points step by step in the following.

A popular model to describe the random variability of quasars is a so-called damped random walk, which is the only process that is stationary, Gaussian, and Markov. This process is also known as an Ornstein–Uhlenbeck (OU) model, which is the nomenclature we will use here. It is similar to red-noise models with a broken power law as the power spectral density, which have been used to model Sgr A* in the past (Do et al. 2009; Meyer et al. 2009; Witzel et al. 2012). Most recently, Dexter et al. (2014) used an OU process to successfully describe the submillimeter variability of Sgr A*. Additionally, Kelly et al. (2009), MacLeod et al. (2010), and Zu et al. (2013) have shown that an OU process is an excellent description of active galactic nucleus (AGN) flux variations. Its key advantage for our purpose is that it is a time continuous model that makes the handling of sampling gaps easier.

The OU model is determined by three parameters, the mean μ, the speed of mean reversion k, and the volatility σ, and the dynamics can be described by the following stochastic differential equation:

$$\begin{equation} dy = k(\mu - y) dt + \sigma dZ_t, \end{equation} \tag{ 9 }$$

where dZ_t is the increment of a Brownian motion with $Z_t \sim \mathcal {N} (0,t)$. A solution to this equation is the conditional distribution

$$\begin{equation} f(y_{t+\Delta t} | y_t) = \frac{1}{\sqrt{2\pi \tilde{\sigma }^2}} \exp \left(\frac{-(y_{t+\Delta t} - \tilde{\mu })^2}{2\tilde{\sigma }^2}\right), \end{equation} \tag{ 10 }$$

$$\begin{equation} \tilde{\mu } = y_t e^{-k\Delta t} + \mu (1- e^{-k\Delta t}), \end{equation} \tag{ 11 }$$

$$\begin{equation} \tilde{\sigma }^2 = \frac{\sigma ^2 (1- e^{-2k\Delta t})}{2k}. \end{equation} \tag{ 12 }$$

Our goal is to identify whether more than one state is present in the light curves from Sgr A*. Past studies of this source have already revealed a great deal of information, and in particular imply that any additional state will be subtle: as reported by Witzel et al. (2012) and Dodds-Eden et al. (2011), a unimodal, heavy-tailed, power-law-like distribution convolved with a Gaussian describes the flux density distribution accurately. In our case, it is therefore important that the baseline model is already a good description of this overall flux distribution and power spectral density; otherwise, a model with two or more states might be wrongly preferred over the baseline model of one state if the one state is modeled incorrectly. An OU process matches the observed characteristics of Sgr A*'s power spectrum: a broken power law that is otherwise featureless. The observed flux distribution, however, is not Gaussian. In order to match the power-law-like distribution convolved with a Gaussian, we will also use the exponential and double-exponential of an OU process to model Sgr A*. Since this translates into the equivalent of a lognormal (and loglognormal) distribution, we will use the notation of a log - and log log -OU process in the following. Please note that taking the logarithm twice is purely empirically motivated and not rooted in physical considerations. Ideally, we would like to use a power-law distribution, but a stochastic model similar to Equation (9) which results in a power-law density is not known.

The conditional distributions for the log - and log log -OU processes are given as follows³: (1) for the log  of a OU process,

$$\begin{equation} f(y_{t+\Delta t} | y_t) = \frac{1}{\sqrt{2\pi \tilde{\sigma }^2}y_{t+\Delta t}} \exp \left(\frac{-(\ln {(y_{t+\Delta t})} - \tilde{\mu })^2}{2\tilde{\sigma }^2}\right), \end{equation} \tag{ 13 }$$

$$\begin{equation} \tilde{\mu } = \ln {(y_t)} e^{-k\Delta t} + \mu (1- e^{-k\Delta t}), \end{equation} \tag{ 14 }$$

$$\begin{equation} \tilde{\sigma }^2 = \frac{\sigma ^2 (1- e^{-2k\Delta t})}{2k}, \end{equation} \tag{ 15 }$$

and (2) for the log log  of a OU process,

$$\begin{eqnarray} f(y_{t+\Delta t} | y_t) & = & \frac{1}{\sqrt{2\pi \tilde{\sigma }^2}y_{t+\Delta t}\ln (y_{t+\Delta t})}\nonumber\\ &&\times \exp \left(\frac{-(\ln (\ln (y_{t+\Delta t})) - \tilde{\mu })^2}{2\tilde{\sigma }^2}\right), \end{eqnarray} \tag{ 16 }$$

$$\begin{equation} \tilde{\mu } = \ln (\ln (y_t)) e^{-k\Delta t} + \mu (1- e^{-k\Delta t}), \end{equation} \tag{ 17 }$$

$$\begin{equation} \tilde{\sigma }^2 = \frac{\sigma ^2 (1- e^{-2k\Delta t})}{2k}. \end{equation} \tag{ 18 }$$

The choice of a time continuous model leads to the parameter Δt in the equations above and therefore offers a direct way to deal with the gaps in Sgr A*'s light curve for the conditional distribution. However, the transition matrix P = p_ij is calculated for a specific time difference τ and must also be modified when gaps are present.⁴ A straightforward way to do this is to multiply the transition matrix with itself N times for a gap that is τ · N:

$$\begin{equation} A = P^N, \ {\rm with }\ N = \max (1, {\rm round}(\Delta t / \tau)). \end{equation} \tag{ 19 }$$

For Sgr A*'s data set, which has an average sampling of one measurement per 1.2 minutes (not counting the large nightly/yearly gaps), we chose a final value of τ = 1 minute. We have explored much shorter values but found it to make no significant difference. For values of Δt > 1000 minutes, we set Δt = 1000 minutes since these gaps are safely greater than the mean-reversion timescale of Sgr A* (Meyer et al. 2009; Witzel et al. 2012), and the flux density points after these gaps are therefore a new, independent realization.

The stochastic model in Equation (9) aims to describe the intrinsic properties of the source under consideration. In a realistic setting, however, an additional noise component is present that reflects the measurement process. This measurement noise is typically white noise, i.e., it does not depend on the previous observation, and it is well described by a Gaussian and therefore fully specified by one parameter σ_meas. Often, there exists an estimate of σ_meas. For the case of Sgr A*, for example, nearby stars of similar magnitude visible in the same image offer a straightforward way to estimate the measurement noise since these stars have constant flux intrinsically. In case an independent estimate of σ_meas is present, it can be advantageous to include it in the stochastic model of the source. In order to do this, the conditional distribution from the OU process has to be convolved with a Normal distribution $\mathcal {N} (0,\sigma _{\rm {meas}})$,

$$\begin{eqnarray} \tilde{f}(y_{t+\Delta t} | y_t;\sigma _{\rm {meas}}) &=& \int f(x|y_t) \cdot \frac{1}{\sqrt{2\pi \sigma _{\rm {meas}}^2}} \nonumber\\ &&\times \exp \left(\frac{-(y_{t+\Delta t} - x)^2}{2\sigma _{\rm {meas}}^2} \right) \; dx. \end{eqnarray} \tag{ 20 }$$

For the important cases of a log - or log log -OU process (Equations (13) and (16)), this integral has to be calculated numerically in every term of the sum in the likelihood function of Equation (8). Note that this dramatically increases the computing time. For us, the duration of computing the posterior with the numerical integration increased to the order of days from just a few hours without it. Please also note that this approach of incorporating the measurement noise assumes that it is constant. While photon noise leads to an increase of noise with flux, the data properties are well modeled with the assumption of constancy and may reflect the dominance of point-spread function measurement noise (Witzel et al. 2012; G. Witzel et al., in preparation).

An important question is how to decide whether more than one state is needed at all, and if so, how many different states can be inferred from the data. Here lies an important distinction between astrophysical and economic analyses: while the latter are mainly interested in a precise model of a time series to make accurate forecasts, the former are insight driven. The necessity of different states could point to different physical mechanisms and elucidate the astrophysics of the source under consideration. Whether or not an additional state is necessary in the model of an astronomical time series should be assessed by several metrics. Bayesian methods such as comparing the Bayesian evidences for different models belong to the standard methods well suited to comparing models but have the (dis)advantage of forcing one to write down specific priors for the parameters, which are often ambiguous. Note that Hamilton (2008) warns that methods relying on likelihood ratio tests fail to satisfy necessary regularity conditions.

A quite robust way of assessing the necessity of an additional state is to look at the global likelihood/posterior distribution in another way: too many assumed distinct states will lead to a highly multi-modal, very irregular looking posterior distribution. If only the maximum of the distribution is determined, the optimizer might often fail to converge at all in that case. Most importantly, the individual states should show persistence, if they are real. In any solution with a superfluous additional state, the probability of remaining in that state p₂₂ will be ≪0.1, while a significant additional state should show persistence with p₂₂ ≳ 0.8. If the additional state has a very low probability of remaining in that state, two things might occur. (1) The probability of remaining in the first state is also very low, meaning that the states fluctuate from point to point. It seems extremely unnatural that state changes in the observed source occur exactly at the sampling rate of the measurements. (2) The probability of remaining in the first state is very high. In this case, the source will be in the additional, second state hardly at all. Only very few flux points will be assigned to that state, and a natural explanation is that these are outliers for the assumed first state conditional distribution. This can easily happen if any assumption of the stochastic model is not quite accurate, e.g., if the flux density distribution is not quite lognormal or the measurement noise is not strictly constant.

Up to now, the timing behavior of Sgr A* has always been modeled using only a single state in the literature. An advantage of this approach is that it allows for any choice of the overall flux distribution, e.g., a power law convolved with a Gaussian as in Witzel et al. (2012). In our methodology, however, where we explicitly model conditional distributions, we do not have the freedom to choose a power-law behavior since a dynamic model similar to Equation (9) that results in power-law distributions is not known. Therefore, we have to show first that a log - or log log -OU process can accurately model the overall, unconditional flux distribution as well. If this was not the case, any additional state would likely be preferred just because a single state does not appropriately reflect the distribution of fluxes.⁵ The key feature we are looking for is a different timing behavior since this is expected from a significantly distinct state that reflects a distinct physical process in the accretion flow.

Most importantly, the total flux distribution of Sgr A* is peaked and highly skewed. Exploring the distribution visually, it is noticeable that the histogram of log (flux) and log (log (flux)) resembles a Gaussian, although it is still skewed in the latter case. This skew is minimized when a constant is added to the entire light curve, i.e., y_t = y_t + c for all t, where y_t denotes the flux density at time t. For the loglog of the flux, this constant is evaluated at c = 1.25 mJy for the VLT data and c = 1.35 mJy for the Keck data, and in the case of log(flux), it is 0 mJy for both data sets. Note that a constant of c > 1 would be required anyway in order for the loglog to be defined since the light curves are normalized such that min (y_t) = 0. Table 1 shows the Bayesian evidences and the best-fit parameter values for both models and data sets. The preferred model in both cases is the loglog one, although the difference from the log model is marginal for the VLT data. We adopt the loglog-OU process as the baseline when measurement noise is not included via Equation (20), and the resulting model as well as the observed flux density histograms for the VLT data are shown in Figure 3. The comparison of baseline model and data shows that while the agreement is very good, it is not perfect. However, since we cannot use power-law distributions and do not explicitly model the measurement noise here, this is not surprising.

Zoom In Zoom Out Reset image size

We have used the Bayesian nested sampler MultiNest (Feroz & Hobson 2008; Feroz et al. 2009) to explore the posterior distribution and assumed uniform priors for the parameters $\mathcal {U}(0, 10)$.

We can now turn to the question we set out to answer: does Sgr A* have more than one state? While the choice of the conditional distribution for the baseline model was motivated to most closely match the observed overall flux distribution (for a single state, the conditional flux distribution obviously equals the overall one), the two-state case is more complex. Fundamentally, the selection of the form of the conditional distributions (e.g., OU, log-OU, etc.) is an assumption to make and cannot be derived from the data. In a practical approach, we have chosen many different combinations and looked at them individually.

The results of the modeling with different assumptions for the conditional distributions are summarized in Table 1. Strikingly, all two-state scenarios are preferred over the single-state model for both data sets as indicated by the Bayesian evidence. All posteriors are well behaved and not more than bimodal (for a two-state model, the posterior will often be bimodal⁶), and the probability of remaining in a given state is high, indicating persistence of these states. Taken together, this is strong evidence for the presence of more than one variability state in the observed flux from Sgr A*. The overall preferred two-state model is one consisting of a log log -OU and a log -OU state. This is true for both the VLT and Keck data sets.

In Figure 4, we show the observed VLT flux density time series in a decomposed way: the upper panel contains all points that are with probability >0.5 in the log -OU state, while the lower panel shows all flux points that are with probability >0.5 in the log log -OU state. The differences are clearly visible. In fact, in all explored two-state models, one state is always reverting to its mean fast and has a lower mean, while the other state is slower mean reverting with a higher mean. A fast mean-reversion here means that this timescale is at the order of the sampling. We have also tested that the second state is not just preferred because of a better description of the overall flux distribution: forcing the two states to have the same form of the unconditional distribution by setting μ₁ = μ₂ and σ₂ = k₂/k₁ σ₁, we still see strong evidence for a second state (see last lines in Table 1). It is the timing that drives the significance of a second state.

Zoom In Zoom Out Reset image size

The fact that one of the two states is quickly mean reverting to a fairly low mean suggests that the interpretation of this state predominantly describes the parts of the light curve that are dominated by measurement noise, which is approximately white noise. Note that even if this is the case, the above analysis is still valuable since it offers a recipe to label each point-to-point change of the flux series as "noise dominated" or "source dominated." Consistent with this interpretation is the observation that the less noisy Keck data are only "noise dominated" 6% of the time, while the VLT data are "noise dominated" 17% of the time.

Our hypothesis is that the quickly mean-reverting state represents the measurement, white noise process. We will test this in two ways. Since we have an estimate of σ_meas, we can make use of Equation (20). Furthermore, Witzel et al. (2012) developed a Monte Carlo tool to simulate Sgr A* light curves (assuming a single state only), and we can analyze these mock data separately.

We summarize our results in Table 2 for the real data and in Table 3 for the mock data. Regarding the real data, the evidence for a second state vanishes when the measurement noise is incorporated into the model. For both data sets, the second state becomes highly non-persistent with a probability of remaining in that state of p₂₂ ≲1%. This is in stark contrast to the case of modeling without the extra noise component, where p₂₂ ≳ 90% generally. For the VLT data, there is also no improvement in the Bayesian evidence. While there is some improvement in the evidence for the Keck data, the structure of the solution suggests that a few outliers drive this behavior (see also the discussion in Section 2.3).

Regarding the mock data, Table 3 shows that the mock data lead to a very similar solution as the real data, which is encouraging. Since the mock data algorithm was calibrated to the VLT data in terms of the sampling function and applied noise, the agreement between the mock data and the real VLT data is stronger than with the Keck data. When the measurement noise is not incorporated into the modeling, a second state is clearly preferred. Since the mock data follow a single-state process by construction, this shows that the white measurement noise leads to a discernible state in addition to the intrinsic red noise process of Sgr A*. Again, when the noise is explicitly convolved into the conditional distributions, the evidence for a second state vanishes. This solution puts the probability of remaining in the second state at p₂₂ ≲ 1%.

In summary, the evidence for two distinct states in the flux from Sgr A* disappears when the measurement noise is taken into account. Therefore, we see our hypothesis verified: the quickly mean-reverting, lower-mean state represents the state when the observed flux changes are dominated by instrumental noise, and the slowly mean-reverting, higher-mean state represents the state when the flux changes are dominated by Sgr A* itself.

It is interesting to note that the same result holds true when the analysis is done on a combined Keck and VLT data set (see Figure 2). While two states are again being picked up when the measurement noise is unaccounted for, one state tends to describe white noise and the other Sgr A* intrinsically. With the measurement noise convolved into the conditional distributions, the evidence for a second state vanishes.

We have shown that Sgr A* is sufficiently described by a single-state process once the measurement noise is accounted for in the modeling. The accurate model for Sgr A*'s light curve is therefore the log-OU process convolved with Gaussian noise (see Table 2). While it is tempting to interpret these parameters and compare them to different modeling approaches of the past, one has to be cautious when the single-state parameter inference is done as the limiting case of a multi-state model, as we have done here.

The reason for caution is as follows. Our approach models the conditional distribution from one measurement to the next, f(y_{t + Δt}|y_t), which is necessary in the context of a multi-state Markov model. If only a single state is assumed, this means, however, that only the typical sampling horizon Δt is used to estimate the mean-reversion timescale k. If the OU process is completely accurate and describes the true nature of the source, this is unproblematic. If, however, there are deviations from the OU process in the data, a different sampling horizon would lead to a different estimate of k. Ideally, a global approach using all time lags in the data would be used to estimate k for a single-state process. Since we focus here on multi-state modeling, this is beyond the scope of this work. We would like to note, however, that the analysis of our mock data lead to consistent results with the analysis of the real data, which means we explicitly confirm the results of Witzel et al. (2012).