香京julia种子在线播放

Except for five roving planets and an occasional comet or nova, the nighttime sky seems immutable to the human eye. The pattern and brightness of stars appears unchanging as from our childhood to old age. Myths from ancient Egyptian, Greek, and Australian Aboriginal cultures suggest that a few stars (such as Algol, Mira, and Aldeberan) were recognized as variables [1, 2]. As telescopic studies proliferated from the seventeenth through twenty first centuries, more variable stars were found with a wide range of characteristics. Some are periodic due to pulsations, rotationally modulated spots, or eclipses of binary companions. Others vary in irregular ways from magnetic flares, eruptions, pulsations, accretion of gas from companions, and most spectacularly, nova and supernova explosions. Ten thousand stars in two dozen categories were cataloged by Kukarkin and Parenago [3]; this catalog now has over 50,000 stars with >100 classes [4]. NASA's Kepler mission has recently shown that most ordinary stars are variable when observed with ~0.001% accuracy and dense cadences [5].

The study of celestial objects with variable brightness has broadened hugely in recent decades, emerging as a recognized discipline called “time domain astronomy” [6]. The brightest sources in the X-ray and gamma-ray sky are highly variable, typically from accretion of gas onto neutron stars and black holes. Timescales range from milliseconds to decades with a bewildering range of periodic, quasi-periodic, stochastic and bursting characteristics. The Galactic black hole binary GRS 1915+105 alone has a dozen modes of variability [7]. The radio sky has extragalactic quasars and blazars as well as Galactic pulsars and several varieties of fast radio bursts and transients. The non-photon gravitational wave observatories have recently emerged with rapid “chirps” from merging black hole and neutron star binaries [8]. A huge industry searching for distant supernova explosions is propelled by their utility in tracing the accelerated expansion of the Universe [9].

Given this tremendous growth in the volume and complexity of astronomical time series data, we can ask what methods are common for their characterization and analysis. In the statistical literature, the dominant methods are time domain parametric modeling. Popular introductory texts [10, 11] devote most of their chapters on autoregressive time domain models with only brief treatment of Fourier frequency domain methods. Prominent graduate level textbooks have a similar balance [12, 13]. Yet only few research papers each year in the astronomical literature mention autoregressive time domain methods: <10 papers/year prior to 2010, recently growing to ~25 papers/yr. In contrast, ~1,000 papers annually use “Fourier” or “power spectrum” procedures.

This radical dichotomy between the methodologies emphasized by statisticians and those used by astronomers motivates our investigation here. Perhaps the difference arises due to real differences in data characteristics and questions raised in a physical science like astronomy compared to questions raised in social sciences like econometrics that are treated in the statistical literature. This seems possible: time series texts oriented toward engineers [14, 15] and meteorologists [16] generally use spectral and wavelet analysis rather autoregressive modeling. Another difference is that social science time series mostly have evenly spaced observation times while astronomical time series are often irregularly spaced due to daily or annual celestial cycles and other causes. The recent growth in autoregressive models for astronomy is mostly restricted to the narrow field of continuous-time models that treat irregular time series [17].

Our goal here is to investigate whether standard time domain methods commonly taught and used in the social sciences can in fact be effective in characterizing astronomical time series, and whether this analysis path can lead to valuable scientific insights.

We start with a brief presentation of the mathematics of common time domain models, model fitting and validation, and implementation of modeling in public domain software. We proceed with application of standard methods to the brightness variations of four stars using both evenly- and irregularly-spaced light curves based on the work of Caceres et al. [submitted for publication, in preparation]. We continue with an outline of various extensions to the common autoregressive models, including new developments of continuous-time models appropriate for irregularly-spaced time series. We end with discussion on what astronomical insights can emerge from autoregressive modeling, and comments on the utility of this approach to future time domain studies in astronomy.

In this section, we give a simplified presentation of standard material on parametric autoregressive modeling. The suite of ARMA-related procedures was brought into prominence by the 1971 textbook of Box and Jenkins [[18] is the latest edition] and is sometimes known as the Box-Jenkins method. Additional coverage appears in time series analysis texts like Chatfield [10] and (in greater mathematical detail) Hamilton [12], and econometric textbooks like Enders [19], Greene [20], and Hyndman and Athanasopoulos [21]. The ARFIMA (also called FARIMA) model is less common than the ARIMA model (~0.05 M vs. ~2 M Google hits). Its mathematical foundations are treated in monographs such as Palma [22]. Autoregressive modeling is reviewed in the astronomical literature by Scargle [23], Koen and Lombard [24], and Caceres et al. [submitted for publication].

As elsewhere in statistical analysis, it is wise to examine the data nonparametrically prior to parametric modeling. The time series should be examined for its basic behavior with kernel density estimation (or similar nonparametric regression) smoothing to see trends. Box-plots, histograms and quantile-quantile plots show the distribution of intensity values. Most importantly here, the nonparametric autocorrelation function (ACF) gives the correlation between a series and time-lagged values of itself over its entire length. For an evenly-spaced time series x with time stamps t, mean value μ and standard deviation σ, (1)ACF(k)=E[(xt-μ)(xt-k-μ)]σ2. In realistic cases, the population mean and standard deviation are unknown and must be estimated from the data. As with Fourier and wavelet transforms, the full information in a time series is reproduced in the ACF if an unlimited number of coefficients are used. Theorems give the statistical distribution of the ACF for Gaussian white noise, so confidence intervals (such as P = 99%) can be constructed and hypothesis tests applied (Durbin-Watson and Breusch-Godfrey tests) to readily evaluate whether autocorrelation is present.

When the ACF shows significant signals, current values of the time series depend on past values. The time series violates the assumption of “independence” in i.i.d. and many standard statistical procedures can not be accurately applied. This is the situation where autoregressive modeling can be effective. Commonly, two forms of dependency on past values are treated as a linear regression. First, an autoregressive (AR) process has coefficients that quantify the dependence of current values on recent past values: (2)xt=a1xt-1+a2xt-2+…+apxt-p+ϵt where ϵ_t is a normally (Gaussian) distributed random error with zero mean and constant variance, p is the order of the process (i.e., how many lags to model), and the a's are the corresponding coefficients for each lag up to order p. Second, a moving average (MA) process has coefficients that quantify the dependence of current values on recent past random shocks to the system: (3)xt=ϵt+b1ϵt-1+b2ϵt-2+…+bqϵt-q where ϵ_t is the error term for the t-th time point, b_i is the coefficient for each lagged error term up to order q. Here the past ϵ noise values can be viewed as random shocks to the system, and the b coefficients quantify the response to the shocks, ϵ. Adding these two equations together gives a combined ARMA(p,q) process. Coefficients are estimated by standard regression procedures such as maximum likelihood estimation.

The ARMA model has several assumptions including stationarity that requires constant mean and constant variance. A stationary time series exhibits the same behavior at the beginning, middle, and end of the dataset. Violations of constant variance (homoscedasticity) can be treated with nonlinear ARMA-type models with volatility like GARCH. Violations of stationarity can be tested and regimes of different behavior can be found using change point analysis.

Nonstationarity and variable mean values can sometimes be removed by fitting a global regression model such as a polynomial, but often an adequate detrending regression model cannot be found. A flexible nonparametric procedure called differencing can remove nonstationarity in many such cases. Here one applies the backshift operator B that replaces the time series x_t by another xt′ consisting of the point-to-point difference in values: (4)xt′ = xt-Bxt = xt - xt-1. The differenced time series can then be modeled as a stationary ARMA process, and the original time series with trends can be recreated by reversing or integrating the differenced time series. This combination of nonparametric differencing and integration with a parametric ARMA process is called the ARIMA(p,d,q) model where d represents the number of differencing operations applied and typically equals one.

Finally, and less intuitively, a fractional integrated procedure can be described by (5)(1-B)dxt=ϵt where d can be a real (non-integer) order of differencing and B is the backshift operator defined above. When combined with ARMA modeling, this gives an ARFIMA(p,d,q) process. The advantage here is that fractional values of d quantifies a long-memory dependency of x_t on past values in a single parameter. Mathematically, this fractional integration process corresponds to the Fourier behaviors known as 1/f^α “red noise” where f is the frequency and α is the slope of a power law fit to the low-frequency Fourier power spectrum [22]. It can also be viewed as a binomial series expansion, (6)(1-B)d=∑k=1∞(nk) (-B)k. The parameter d is typically constrained to be in the range 0.0 − 0.5 to define a stationary model. Three commonly used parameters are related to each other as follows: d = α/2 and d = H − 1/2 where H is the econometricians' Hurst parameter.

In simple language, ARFIMA(p,d,q) is a time series model where the autoregressive and moving average components treat many short-memory dependencies on recent past values, and the fractional integration component treats both many forms of trend and a 1/f^α-type long-memory dependency. Best fit parameters for ARMA -type coefficients are typically calculated with maximum likelihood estimation for different values of p, d, and q, called the order of the model. The optimal order is then selected to optimize the Akaike Information Criterion (AIC), a penalized maximum likelihood measure that balances improvements in likelihood with increases in model complexity.

ARMA, ARIMA, and ARFIMA models are attractive for astronomical time series analysis for various reasons. First, the models are very flexible, successfully modeling an astonishing variety of irregular or quasi-periodic, smooth or choppy, constant or variable mean light curves. Second, the dimensionality of the model is relatively low with a moderate computational burden of the numerical optimization. In contrast, local regression methods, such as Gaussian Processing regression, may have thousands of parameters with O(N³) computational burden although some computational efficiency has been achieved [25]. Third, error analysis on the parameters naturally emerges through the likelihood regression analysis. Fourth, they are extensible to situations involving multivariate time series, combinations of stochastic and deterministic behaviors, change points, and (moderately) irregular observation spacing. Autoregressive modeling is not well-adapted to situations with strictly periodic variations (where the signal is compactly concentrated in Fourier power coefficients) or with sudden eruptive events (where the nonstationary amplitude is not greatly reduced by differencing).

Altogether, ARIMA-type models have proved to be useful to a very wide range of physical and human-generated dynamic systems in engineering, econometrics, Earth sciences, and other fields. But despite these advantages, non-trivial ARIMA models have been used only rarely in time domain astronomy for compact radio sources [26], Mira stars [27] variability, and quasars [28]. ARFIMA models appear in only a single study on solar flare activity [29].

The first action that every analyst should take after ARIMA-type modeling is model validation based, in part, on residual analysis. This should follow any regression procedure for two reasons. First, even though the “best fit” has been obtained in a maximum likelihood sense, the entire model family may not apply to the dataset under study; this problem is called “model misspecification.” Second, the model may be correctly specified but its underlying mathematical assumptions may be violated. For example, many regression procedures assume the errors (ϵ in Equation 1) are “normal i.i.d.”; that is, they are independently and identically distributed following a Gaussian distribution.

As with any regression analysis, the unsequenced values of model residuals can be tested for normality using the nonparametric Anderson-Darling or Shapiro-Wilks tests. The Jarque-Bera test for normality examined the skewness and kurtosis of the residual distribution. The Breusch-Pagen and White tests are used to test whether the residuals have constant variance or whether heteroscedasticity is present; this is a violation of the “identically” requirement in i.i.d.

A suite of specialized “time series diagnostics” have been developed by econometricians that use the time sequence of residuals to test model validity in various ways. These are hypothesis tests for specific characteristics where the statistic has a known distribution for the null hypothesis of Gaussian white noise. Thus quantitative probabilities for deviations from Gaussian white noise can be obtained under many circumstances.

The most well-known of these diagnostics is the Durbin-Watson test of serial autocorrelation; i.e., whether the AR(1) coefficient of the residuals is consistent with zero. The Breusch-Godfrey test generalizes this to test significance of AR(p) where p is specified by the user. The Box-Pierce and Ljung-Box test are “portmanteau” tests that address all values of p simultaneously. Other nonparametric tests (such as the Wald-Wolfowitz runs test that groups the residuals into a binary “+” and “−” variable) may also be effective in testing the existence of temporal structure in model residuals.

Several statistics test for nonstationarity in the residuals. A simple nonparametric correlation coefficient like Kendall's τ can test for overall correlation. Global regressions can reveal large-scale trends and local regressions smooth the residuals and reveal small-scale structure. But it is possible for a time series to be nonstationary but without trend. One of these conditions is called “unit root” where the time series does not revert to the original mean after experiencing a shock. The physicists' random walk is an example of a nonstationary diverging process. The augmented Dickey-Fuller and Kwiatkowski-Phillips-Schmidt-Shin (KPSS) statistics test for trend stationarity against a unit root.

Software for autoregressive modeling is widespread. In Python, statsmodel.tsa and related libraries have ~70 functions implementing AR, ARMA, ARIMA, and VAR (vector autoregressive) model fitting. The kalmanf function estimates ARMA models with exact maximum likelihood estimation using the Kalman filter. The regime_switching function implements a nonlinear Markov switching dynamic regression. Time series diagnostics, trend regressors and other ancillary functions are included.

Matlab has several dozen infrastructure functions for manipulating time series objects. But its power is in the Econometrics Toolbox with extensive methodology (and tutorials) including: time series decomposition, filtering, ARIMA regression with non-i.i.d. errors, ARMAX estimation with exogenous covariates, GARCH models, multivariate VAR, AIC model selection, state space modeling with the Kalman filter, time series diagnostics, simulation and forecasting. This Matlab software costs several tens of dollars for a student, several hundred dollars for a university researchers, and several thousand dollars for a commercial license.

The largest collection of autoregressive software resides in the public domain R statistical software environment [30]. Base-R software includes infrastructure for “ts” (evenly spaced) and “zoo” (irregularly spaced) time series objects as well as basic ARIMA and ARIMAX functionality. But the bulk of advanced time series methods are implemented in some of the >12,000 add-on CRAN packages. Mature and extensive CRAN packages for ARIMA-type modeling include acp, arfima, BAYSTAR, BigVAR, bsts, carfima, carx, cents, cts, dse, fGarch, FGN, FitAR, FitARMA, forecast, glarma, imputeTS, ltsa, MARSS, mclcar, mlVAR, MSBVAR, MTS, partsm, perARMA, portes, prophet, robustarima, rugarch, sparsevar, stochvol, stsm, TSA, tsdecomp, tsDyn, tseries, and vars. Recently, the continuous-time CARFIMA model has been implemented in CRAN package carfima. Together these packages provide more than a thousand functions relating to autoregressive modeling. The package capabilities are briefly listed and classified in the CRAN Task View on Time Series Analysis. An elementary introduction to R for astronomical time series appears in the text by Feigelson and Babu [31].

Important diagnostics for goodness-of-fit for time series and model fit residuals (such as the Durbin-Watson, Bruesch-Godfrey, Box-Pierce, Augmented Dickey-Fuller, and KPSS tests) are available in CRAN packages such as tseries [32] and lmtest [33]. These tests are described in econometrics textbooks [e.g., [19, 20]] and in Wikipedia.

The ARIMA and ARFIMA fits in Figures 1–4 here were obtained using the auto.arima and arfima functions in forecast, one of the most popular CRAN packages that is describe in the ebook by Hyndman and Athanasopoulos [21]. The data manipulation and plotting was performed using base-R. Benchmarking shows that the auto.arima and arfima codes provide AIC-optimized maximum likelihood fits in O(N) wall-time: ~0.4 CPU-sec for ≤1,000 datapoints, in 4 CPU-sec for 10,000 datapoints, and so forth. This is sufficiently efficient for analyzing large astronomical surveys. However, the new carfima code (below) is much slower.

In addition to ARFIMA for 1/f^α noise processes, there are many mathematical variants of ARMA and ARIMA modeling. They are typically developed for economics time series analysis and are described in econometric textbooks like Enders [19], Greene [20], and Hyndman and Athanasopoulos [21].

A mathematical limitation of classical autoregressive modeling is the restriction to evenly spaced datasets. However, astronomers need methods that are robust to irregularly spaced observing cadences limited by daily and annual celestial cycles as well as telescope allocation constraints. The underlying signal is typically generated by a continuous process of complex, often nonlinear, physical mechanisms that are not fully understood such as accretion disks and magnetohydrodynamical turbulence. Analysis is commonly based on Fourier analysis, but its underlying theorems apply only for the limited situation of a sinusoidal signal superposed on Gaussian white noise. This model is often inadequate to characterize the often multi-scale nonstationary structure in the time series. Finally, a variety of analysis techniques are needed to capture different properties in order to classify heterogeneous ensembles of time series, such as the 40 billion unevenly spaced light curves expected from the forthcoming Large Synoptic Survey Telescope.

In the stellar light curve applications above, we show that a sufficiently dense cadence of irregular observations can be successfully binned or interpolated onto a fixed time grid with “missing data,” allowing classical ARIMA-type modeling to proceed. However, a more effective approach may be generalizations of stochastic difference equations to stochastic differential equations. These have been studied mathematically starting in the 1940s to model physical phenomena such as Brownian motion and econometric variables [41, 42]. Just as ARMA processes play a central role in the representation of time series with evenly spaced observations, CARMA (continuous-time autoregressive moving average) processes play an analogous role in the representation of time series with continuous time parameter allowing direct treatment of irregular time sampling. Many variants are studied: Brownian motion and fractional Brownian motion with long-range dependence, Ornstein-Uhlenbeck process with exponential damping, Lévy process with heavy-tailed and asymmetric time series. The interest in CARMA-type processes is motivated mainly by the successful application of stochastic differential equation models to problems in finance.

In astronomy, a continuous autoregressive CAR(p) process was first applied by Koen [43]. A simple CAR(1) model, variously known as the Ornstein-Uhlenbeck process and damped random walk model, was brought into use for modeling quasar light curves in irregular astronomical time series by Kelly et al. [17] and generalized to CARMA(p, q) by Kelly et al. [28]. The latter paper gives the power spectral density, and autocovariance function at lag τ for the CARMA model. Kelly's software implementation called carma_pack on Github provides a Markov chain Monte Carlo sampler, an adaptive Metropolis algorithm combined with parallel tempering, for performing Bayesian inference on these continuous time autoregressive moving average models. This model has been used in dozens of astronomical studies, particularly for study of variability of quasars for which unevenly spaced observations are commonly available [e.g., [44, 45]].

CARMA models have favorable mathematical properties for astronomical use, although the statistical theory of continuous stochastic processes can be difficult [46]. Mathematically, the power spectrum density of a CARMA process is a sum of Lorentzian functions with centroid, widths, and normalizations of the Lorentzian functions as free parameters [47]. This provides a significant amount of flexibility for modeling complicated power spectra, such as accretion disk systems with quasi-periodic oscillations and red noise, suggesting that CARMA modeling may be applicable to many classes of astronomical variables. Moreover, CARMA models have a state space representation enabling the use of the Kalman filter for calculating their likelihood function. The computational complexity of calculating the likelihood function for CARMA models can therefore scale linearly with the number of data points in the light curve, an advantage for application to massive time-domain data sets.

The CARMA models account for irregular sampling and measurement errors, and are thus valuable for quantifying variability, forecasting and interpolating light curves, and variability-based classification. But they suffer the same limitations as ARMA(p, q) for regularly spaced data: it assumes (weak) stationarity of behavior across the light curve (for instance, no trend in brightness) and only treats short-memory autocorrelation out to max(p, q) time intervals. ARFIMA, on the other hand, treats nonstationarity and a powerlaw model of long-memory behavior in addition to these elements.

The mathematical development of continuous-time generalizations of ARIMA and ARFIMA have only recently begun; the effort is reviewed by Brockwell [48]. Theorems underlying CARIMA are presented by McElroy [49] and for CARFIMA (also called FICARMA) by Brockwell and Marquardt [50] and Brockwell [48]. A maximum likelihood estimation procedure for CARFIMA is developed by Tsai and Chan [51] with R public software implementation released by Tak and Tsai [52]. We encourage time domain astronomers to try modeling with their carfima R/CRAN package when faced with irregular light curves of variability that includes trend and long-memory “red noise” characteristics. It will be particularly valuable to compare the performance of standard ARIMA with “missing data” (above) and continuous-time CARIMA on real and simulated astronomical light curves.

Most often, the goal of the econometrician is to forecast future values of the time series. But the astronomer is rarely concerned about forecasting brightness of their target star or quasar. (An exception is solar activity that affects our terrestrial environment.) The astronomer has other goals:

Astronomers have tended to use a narrow suite of time series methods, particularly Fourier analysis. This is appropriate when the physical phenomenon is strictly periodic with roughly sinusoidal shape such as radial velocity of circular Keplerian orbits or acoustic oscillations. The power spectrum concentrates periodic signals but not stochastic behaviors such as autoregression. There is real danger that inaccurate or even incorrect inferences will be made when autoregressive behaviors are present but are ignored by simplistic methods that assume the interesting signal is embedded only in Gaussian white noise. For example, detection of sudden flare events with a “3-sigma” criterion typically assumes Gaussian white noise whereas the effective sample size is reduced for autocorrelated signals. The Bayesian Blocks procedure [62] for detecting intensity changes in count data similarly assumes Poisson white noise.

Over several decades, researchers in many fields have found that low-dimensional autoregressive ARIMA-type models that quantify dependence on past behavior are amazingly effective in modeling many forms of variability, far more flexibly than low-dimensional (e.g., global polynomial) models that involve dependence on time alone. Extensive and user-friendly software is available for ARIMA-type modeling, particularly in the public domain R and Matlab statistical software environments.

Parametric modeling with any mathematical formalism is only useful for time series if the variability behavior is encompassed by the chosen mathematical model. In statistical parlance, the model must be correctly specified. But this constraint can be checked. Once a mathematical family of models is chosen and fitted to the data by maximum likelihood estimation, diagnostic tests can be applied to confirm that the model is indeed successful; for instance, testing that the model residuals are simple Gaussian white noise. Caceres et al. [in preparation] find that the aperiodic variability seen in most stars observed by NASA's Kepler mission can be effectively treated by ARIMA-type models.

A particular limitation of autoregressive modeling is the insistence on evenly spaced observation times. However, in our examination of stellar light curves above, we converted time series with irregular observation times into a sequence of regular intervals by simple binning with “missing data” in empty time slots. Density estimation (smoothing) or imputation procedures can also be applied; ARIMA modeling itself can be used for imputation, filling in missing data in astronomical time series. An alternative approach is to use continuous-time autoregressive models such as CARMA that treats irregular observation times directly with continuous time models.

Our emphasis here is on situations where a significant component of the variation exhibit autoregressive behaviors so that future values depend on past levels and changes. These components could arise as nuisance effects such as atmospheric fluctuations at the telescope, or they could be intrinsic to the astrophysical systems under study. While our treatment has concentrated on stellar photometric light curves, the methodology can readily be applied to other situations arising in astronomical research: optical or X-ray brightness variations from accreting supermassive black holes (quasars and Seyfert galaxies); complex signals in gravitational wave detectors; rapid atmospheric variations in ground-based active optics instruments; space weather induced background in satellite observatories; and so forth.

It is important to recognize that the scope of modern time series methodology is vast, propelled by complex problems in econometrics, signal processing, engineering, and Earth science systems. We have covered only a small portion of the field here. For example, change point analysis could help in study of astronomical transients, noise can be removed from astronomical signals with wavelet or sparsity methods, and irregularly sampled astronomical time series can be smoothed with local regression procedures. We encourage astronomers to be more adventuresome in investigating the utility of sophisticated time series methods, including ARIMA-type autoregressive modeling, to address problems in the burgeoning field of time domain astronomy.

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.