taufit: model and extract the damping timescale in AGN light curves using celerite¶

Walkthrough demonstrating how to simulate and fit a damped random walk (DRW) light curve, using the DRW model kernel in celerite:

$$k(t_{nm}) = a \exp{(-c\ t_{nm})}$$ $$= 2\sigma^2 \exp{(-t_{nm}/\tau_{\rm DRW})}$$

which corresponds to the structure function:

$$SF^2 = 2\sigma^2(1-e^{-t_{nm}/\tau_{\rm DRW}})$$

where $\tau_{\rm DRW} = \frac{1}{c}$ and $\sigma=\sqrt{a/2}$.

Import package methods to simulate and fit DRW

Generate and fit simulated DRW light curve¶

Define function to simulate DRW light curves with similar amplitude fluctuations at a given timescale, typical mean magnitude, and error bars, with Gaussian white noise.

Simulate $\tau_{\rm{DRW}}=100$ days:

Shown above: The power spectrum (top left), light curve and posterior model prediction $1\sigma$ error ellipses (top right), $\tau_{\rm{DRW}}$ posterior (bottom left) with 16th and 84th percentiles marked as dashed lines, and the autocorrelation function of the squared residuals (ACF $\chi^2$; bottom right). The ACF $\chi^2$ should located within the grey region (consistent with white noise).

Important: The grey shaded region in the $\tau_{\rm{DRW}}$ posterior plot (bottom left panel) corresponds to > 20% of the light curve baseline. If the distribution peaks in this region, the light curve is not long enough and the result is unreliable (see Kozlowski 2017).

Use corner to plot the posteriors with celerite's native variables or the equilivlant structure-function-convention variable:

Celerite has been shown to be ${\sim}10 \times$ faster than codes such as carma_pack, while having equivilent model PSDs to CARMA (Foreman-Mackey et al. 2017).

Prediction¶

Posterior prediction at future times:

Simulate complex kernel (CARMA model)¶

The kernel function for a general CARMA$(p, q)$ model is

$$k(t_{nm}) = \sum^{p}_{j=1} A_j \exp(r_j t_{nm}).$$

In celerite, this corresponds to a summation of ComplexTerms:

$$= \sum^{J}_{j=1} \frac{1}{2}\left[(a_j+ib_j)e^{-(c_j+id_j) t_{nm}} + (a_j-ib_j)e^{-(c_j-id_j) t_{nm}} \right]$$

where each CARMA term corresponds to a celerite term with $a_j=2\ \text{Re}(A_j)$, $b_j=2\ \text{Im}(A_j)$, $c_j=-\text{Re}(r_j)$, $d_j=-\text{Im}(r_j)$ (Foreman-Mackey et al. 2017; Footnote 17). Here, CARMA(p,q) the order is $p=J$ and $q=p-1$. Note, there is no way to explicitly constrain $q$ to be less than $p-1$ using the celerite solver, so we must adopt this general form.

Summing to $J$, we can easily see that if $p=1, q=0$, we recover the DRW kernel (after noting that the imaginary terms $b_j=d_j=0$):

$$= \sum^{1}_{j=1} \frac{1}{2}\left[(a_j+ib_j)e^{-(c_j+id_j) t_{nm}} + (a_j-ib_j)e^{-(c_j-id_j) t_{nm}} \right]$$ $$= \frac{1}{2}\left[ a_1e^{-c_1 t_{nm}} + a_1e^{-c_1 t_{nm}} \right]$$ $$= a_1e^{-c_1 t_{nm}}$$

Note there is a typo in the published version of the paper where the summations are shown incorrectly to sum to 2J.

We recommend you set the bounds and init arguments yourself for better results. Note that the parameters of higher order CARMA-like models quickly become difficult to physically interpret.

Equivilance to CARMA parameters (coming soon?)

Assess bias in the fitting¶

Generate random DRW light curves with a range of input $\tau_{DRW}$, both seasonally-gapped and not, and see how well we can recover the input for each.

You may togggle whether to plot the output or not by changing the default arugments in simulate_and_fit_drw above.

Plot the results¶

Plot input $\tau_{\rm{DRW}}$ vs. output $\tau_{\rm{DRW}}$. Note the confirmation of the bias reported by Kozlowski (2017) when the input $\tau_{\rm{DRW}}$ approaches 10% of the light curve baseline.