PaperlayerView exact revision on arXiv

Source-linked AI summary

Flexible and Scalable Methods for Quantifying Stochastic Variability in the Era of Massive Time-Domain Astronomical Data Sets

Brandon C. Kelly, Andrew C. Becker, Malgosia Sobolewska, Aneta Siemiginowska, Phil Uttley

arXiv:1402.5978v2astro-ph.IM

TL;DR

Massive astronomical time-domain surveys make quantitative characterization of variability a central analysis problem. This paper presents CARMA modeling for irregularly sampled light curves with measurement errors and finds that it generally describes the examined sources well, except for the RR-Lyrae light curve.

  • Problem

    The central analysis problem in massive time-domain astronomy is quantitatively characterizing variability in the large datasets produced by surveys.

    Source: §1 · p. 1

  • Method

    The paper presents CARMA models for estimating light-curve power spectra, with a likelihood calculation that scales linearly with the number of data points.

    Source: §6 · p. 16

  • Results

    CARMA models generally described the examined light curves well, with the RR-Lyrae light curve as the exception.

    Source: §6 · p. 17

  • Takeaways & Limitations

    CARMA models may provide power-spectrum estimates and variability features for selection and potentially for identifying new classes of variables.

    Source: §6 · p. 16

    Complex, multimodal likelihood spaces can make optimization and sampling problematic, especially for higher-order models and light curves with regular variations.

    Source: §6 · p. 17

Abstract

from arXiv · show

We present the use of continuous-time autoregressive moving average (CARMA) models as a method for estimating the variability features of a light curve, and in particular its power spectral density (PSD). CARMA models fully account for irregular sampling and measurement errors, making them valuable for quantifying variability, forecasting and interpolating light curves, and for variability-based classification. We show that the PSD of a CARMA model can be expressed as a sum of Lorentzian functions, which makes them extremely flexible and able to model a broad range of PSDs. We present the likelihood function for light curves sampled from CARMA processes, placing them on a statistically rigorous foundation, and we present a Bayesian method to infer the probability distribution of the PSD given the measured lightcurve. Because calculation of the likelihood function scales linearly with the number of data points, CARMA modeling scales to current and future massive time-domain data sets. We conclude by applying our CARMA modeling approach to light curves for an X-ray binary, two AGN, a long-period variable star, and an RR-Lyrae star, in order to illustrate their use, applicability, and interpretation.

1. INTRODUCTION

Massive time-domain surveys make quantitative characterization of variability a central challenge, especially for irregularly sampled stochastic light curves. The paper presents CARMA models as a flexible way to model their power spectra while retaining likelihood computation that scales linearly with data points.

1. INTRODUCTION

  • Time-domain surveys will produce millions to hundreds of billions of photometric data points, with LSST expected to collect ∼30 trillion measurements on ∼40 billion objects.These datasets enable systematic multi-passband variability studies and may support discovery of new variability classes and anomalous outliers.Source: §1 · p. 1
  • Quasi-periodic and stochastic light curves are commonly characterized by power spectra or structure functions, with power laws and Lorentzian sums used as parametric forms.The power spectral density describes variability power per frequency interval, while the structure function describes variability amplitude by timescale.Source: §1 · p. 1
  • Irregular sampling can distort periodogram and empirical structure-function estimates, risking false astrophysical differences and variability misclassification.Such distortions can arise because the observed light curve samples a continuous-time stochastic process.Source: §1 · p. 2
  • Monte Carlo forward modeling of periodograms is flexible but computationally expensive, especially for complex models and massive datasets.Uneven intervals of fine and sparse sampling can require very dense simulations or separate periodogram calculations for segments.Source: §1 · p. 2
  • Time-domain Gaussian-process likelihoods account for irregular sampling and measurement errors, but covariance-matrix inversion generally costs O(n^3).This computational cost can be prohibitive for large samples of even moderately dense light curves.Source: §1 · p. 2
  • CARMA models represent power spectra as flexible sums of Lorentzians and retain O(n) likelihood computation through a state-space Kalman-filter representation.The Lorentzian centroids, widths, and normalizations provide flexibility beyond the simple zero-centered Lorentzian of CAR(1).Source: §1 · pp. 2–3

2. CONTINUOUS-TIME AUTOREGRESSIVE

The section defines CARMA(p,q) processes and explains how their roots and coefficients shape their power spectra and autocovariances. It also presents a state-space representation that supports efficient likelihood calculation.

2.1. Definition, Power Spectrum, and Autocovariance Function

  • Stationary CARMA(p,q) processes require q < p and autoregressive roots with negative real parts; minimum phase additionally requires moving-average roots with non-positive real parts.For minimum-phase processes, the input white noise can be uniquely determined from the output process.Source: §2.1 · p. 3
  • CARMA(p,q) autocorrelations are weighted sums of exponentials, yielding exponentially damped sinusoids for complex roots and exponential decays for real roots.Decay timescales and oscillation frequencies are set by the real and imaginary parts of the autoregressive roots.Source: §2.1 · pp. 3–4
  • CARMA(p,q) power spectra are weighted sums of Lorentzian functions whose widths and centroids depend on autoregressive roots and whose weights depend on moving-average coefficients.Lorentzian widths are proportional to |Re(r_k)| and centroids to |Im(r_k)/2π|.Source: §2.1 · p. 4

2.2. State Space Representation

  • The state-space representation describes CARMA dynamics with a p-element latent state vector containing a process and its derivatives, enabling efficient likelihood calculation for measured time series.The solution separates deterministic evolution from the initial condition and a stochastic contribution.Source: §2.2 · p. 4

3. STATISTICAL INFERENCE FOR CARMA

The section develops likelihood-based and Bayesian inference for Gaussian CARMA models, using Kalman recursions to handle sampled time series efficiently. It also describes interpolation, fit diagnostics, model-order selection, and computational strategies for posterior sampling.

3. STATISTICAL INFERENCE FOR CARMA

  • For Gaussian CARMA processes, the likelihood supports maximum-likelihood estimation and Bayesian inference, with autocovariance and measurement-error variances entering the model.The authors assume uncorrelated measurement errors in the time series.Source: §3 · p. 4

3.1. Likelihood Calculation via the Kalman Recursions

  • Kalman filtering reduces CARMA likelihood computation from O(n^3) covariance inversion to O(n) operations for n observations.The method discretizes the continuous-time state-space model at the observation times and incorporates measurement error; the likelihood calculation assumes uncorrelated measurement errors.Source: §3 · p. 4, §3.1 · p. 5

3.2. Interpolation and Extrapolation from a Measured Time Series

  • Conditional Gaussian-process distributions and Kalman filtering provide interpolated or extrapolated values and their uncertainties from a measured time series.For forecasting and backcasting, the specified conditional formulas omit the terms used for interpolation.Source: §3.2 · p. 6

3.3. Bayesian Inference

  • Bayesian inference quantifies uncertainty in CARMA parameters and the inferred PSD, addressing likelihoods that can have multiple maxima for p > 1.The model constrains measurement-error scaling to 1/2 < ν < 2, reflecting an assumed 10% uncertainty in the overall normalization of measurement-error standard deviations.Source: §3.3 · pp. 6–7

3.4. Assessing the Fit

  • Gaussian CARMA fit assessment checks whether standardized Kalman residuals are normal white noise and whether their autocorrelations reveal remaining structure.Autocorrelation in squared residuals can indicate time-varying variance; residuals are evaluated using their sequence index as a regular grid.Source: §3.4 · p. 7

3.5. Choosing the Order of the CARMA Model

  • The authors choose CARMA order by minimizing AICc, since CARMA models are generally non-nested and order selection is framed as model selection rather than hypothesis testing.For a scientifically meaningful PSD feature such as a possible QPO, its posterior quality-factor distribution can assess significance even when including it improves AICc.Source: §3.5 · pp. 7–8

3.6. Computational Considerations and Software

  • Although likelihood evaluation is O(n), its Kalman-filter calculation is serial; multimodal posteriors for p > 1 motivate parallel tempering for MCMC sampling.The method uses hotter chains to move between modes and swaps their parameter values with cooler chains; the authors generally found K = 10 chains sufficient.Source: §3.6 · p. 8

4. EXAMPLE APPLICATIONS TO SIMULATED

Simulated light curves illustrate CARMA modeling under regular and irregular sampling, as well as the behavior of a stationary fit to a process that changes over time.

4.1. Stationary Process Under Regular Sampling

  • 5.04 ± 0.08 days was the estimated timescale of the weak QPO recovered in the selected CARMA(5,1) PSD.AICc selected CARMA(5,1) although the simulated process was CARMA(5,3), and standardized residuals were consistent with Gaussian white noise. Tight PSD constraints below the measurement-noise level reflect extrapolation from the assumed model form, not actual uncertainty there.Source: Figure 2 · p. 9, §4.1 · p. 9, Figure 3 · p. 10

4.2. Stationary Process Under Irregular Sampling

  • The fitted CARMA(5,1) also provided interpolated and forecast light-curve values with 1σ uncertainties, enabling conditional simulation at unobserved times.The simulated irregular light curve used three 90-epoch observing seasons separated by 180 days, with observation spacings of 1–3 days.Source: §4.2 · pp. 9–11, Figure 4 · p. 10
  • With 20% measurement errors and irregular seasonal sampling, CARMA(5,1) captured PSD variability above noise, but the high-frequency QPO below noise was not detectable.The model had the minimum AICc; the authors recommend using the simplest model adequate for variability features above the measurement-noise level.Source: §4.2 · pp. 9–10, Figure 5 · p. 10

4.3. Non-Stationary Process Under Irregular Sampling

  • A stationary CARMA(5,2) fit to the switched process was not flagged by residual diagnostics and produced a PSD blending the two regimes, emphasizing high-signal-to-noise features.The blended PSD picked up the strong QPO in the first half and long-timescale broad-band power in the second; inferred contributions are weighted by signal-to-noise across the sampled frequency range.Source: Figure 8 · p. 11, §4.3 · p. 11, Figure 7 · p. 11

5. EXAMPLE APPLICATIONS TO REAL

The paper applies CARMA models to X-ray binary, AGN, and variable-star light curves, illustrating PSD inference, interpolation, and model limitations across different sources. The examples include densely sampled and highly irregular data, as well as regularly periodic stellar variability.

5.1. X-ray Binary

  • CARMA(5,4) recovers the dominant X-ray-binary QPO at 2.84±0.01 Hz (Q = 16.32 ± 3.27) but misses the weaker ≈6 Hz QPO near the down-sampled frequency limit.AICc selected p = 4, q = 3, which does not include enough QPO features; the CARMA(5,4) residuals were consistent with white noise.Source: §5.1 · p. 12, Figure 9 · p. 12

5.2.1. Kepler Lightcurve

  • The Kepler Zw 229-15 CARMA(6,4) PSD tracks the periodogram and is approximately 1/f^3 below the measurement-noise level, with evidence for flattening to ∼1/f^2 on timescales ≳10 days.The model's standardized residuals showed no evidence of deviations from Gaussian white noise; the PSD resembles a power law with slope ∼−3 on timescales shorter than ≈1 month.Source: §5.2.1 · pp. 12–13, Figure 12 · p. 13

5.2.2. SDSS Stripe 82 Lightcurve

  • For the irregularly sampled Stripe 82 quasar, the CARMA(4,1) residuals are consistent with white noise, while the uncertain PSD is consistent with a power law 1/f^α, α = 2–4.The estimated PSD shows some evidence of steepening toward higher frequencies, but uncertainties are large.Source: §5.2.2 · p. 13, Figure 14 · p. 14

5.3.1. Long-Period Variable

  • For the long-period variable, CARMA(6,0) estimates a characteristic timescale τ = 187 days (95% credibility interval 15.5 < τ < 320 days) and a significant QPO B at 15.9 days.QPO A, at approximately 25 days if real, has only 63% posterior probability and is not statistically significant; QPO B has posterior median quality factor Q = 6.85 (95% credibility interval 1.5 < Q < 42).Source: §5.3.1 · pp. 14–15, Figure 16 · p. 15

5.3.2. RR-Lyrae

  • For the RR-Lyrae, CARMA(7,0) identifies significant mode A at 2.49 days and mutually exclusive modes B or C at 1.30 or 0.56 days, respectively.Mode A is present in 99.986% of MCMC samples; B or C is present with 98% probability, with C more likely. Although residual autocorrelations are consistent with white noise, their distribution is narrower than normal, and convergence was problematic, with different runs yielding different results.Source: §5.3.2 · pp. 15–16, Figure 18 · p. 16

6. DISCUSSION

The discussion presents CARMA models as flexible, statistically rigorous tools for estimating variability in irregularly sampled astronomical light curves, while noting challenges in inference and potential extensions. Applications generally fit the examined sources, with RR-Lyrae light curves as the exception.

6. DISCUSSION

  • CARMA models handle irregular sampling and measurement errors, represent PSDs as sums of Lorentzians, and enable rigorous likelihood-based inference with computation scaling linearly in light-curve length.The paper identifies flexible power-spectrum estimation and variability features for source selection or potentially identifying new variable classes as primary uses.Source: §6 · pp. 16–17
  • Prior extensions include CAR(1) mixtures as CARMA special cases and alternatives to linear Gaussian CAR(1), with Andrae et al. (2013) finding Gaussian CAR(1) best for most Stripe 82 quasars.Kelly et al. (2011) showed their mixtures can approximate broken-power-law PSDs when slopes lie in [−2, 0] and relate to solutions of the stochastic linear diffusion equation; Andrae et al. (2013) also explored non-Gaussian and nonlinear models, while noting CARMA(1,1) is not stationary.Source: §6 · p. 16
  • Complex, multimodal likelihoods—especially for higher-order models and regular variations—are CARMA’s primary shortcoming, requiring care to locate and sample dominant posterior modes.Parallel tempering may still fail with too few chains or insufficient run time; the paper’s maximum-likelihood method uses 100 random starts and local greedy gradient-based optimization, motivating more effective inference algorithms.Source: §6 · pp. 16–17
  • CARMA models generally described the X-ray binary, two AGN, and long-period variable well; for the RR-Lyrae star, they identified the catalogue period in ≈75% of MCMC samples despite the poor fit.The RR-Lyrae residuals deviated from a normal distribution, suggesting residual distributions may help discriminate variable types and that CARMA-derived variability features may still be useful for non-stochastic light curves.Source: §6 · p. 17
  • Potential extensions include modeling periodic light curves with a deterministic component and CARMA residuals, and using multivariate CARMA to characterize multi-passband variability.The paper connects the periodic approach to PARMA models and identifies multivariate CARMA as future work.Source: §6 · p. 17

APPENDIX KALMAN FILTER FOR A CARMA PROCESS

The appendix gives a rotated-state Kalman-filter algorithm for computing conditional means and variances for CARMA time series, then describes how to recursively obtain coefficients for interpolation and extrapolation. These filter outputs also provide the quantities needed to calculate the likelihood.

APPENDIX KALMAN FILTER FOR A CARMA PROCESS

  • The Kalman filter computes each observation’s conditional mean and variance from earlier measurements by initializing the rotated state at its stationary mean and covariance, then iteratively predicting and updating it.The series is centered before filtering, and the first measurement uses the stationary CARMA values.Source: APPENDIX KALMAN FILTER FOR A CARMA… · pp. 18–19
  • The resulting conditional means and variances can be used to efficiently calculate the CARMA likelihood, with the original-series conditional mean recovered by adding μ.The centered-series conditional mean is shifted by μ to obtain the mean conditional on the original observations.Source: APPENDIX KALMAN FILTER FOR A CARMA… · p. 19

ALGORITHM FOR COMPUTING THE COEFFICIENTS NEEDED FOR INTERPOLATION AND EXTRAPOLATION FROM A MEASURED TIME SERIES UNDER THE CARMA MODEL

  • Interpolation and extrapolation coefficients are computed recursively after filtering to the relevant time, with Kalman gains determined by observation times rather than measured values.The algorithm initializes coefficients at t0 and updates them across subsequent observations; if t0 is beyond the final observation, no further calculation is needed.Source: APPENDIX KALMAN FILTER FOR A CARMA… · pp. 19–20
Loading 1402.5978v2…