1. Introduction

Recently, the star Betelgeuse was in the news for sudden dimming, raising excitement that it will go supernova.
Contrary to common belief, the brightness of stars is not constant but can vary over the course of months, weeks, or even hours. Fluctuations in the stellar atmosphere, convection currents, stellar winds, etc. cause both systematic and random changes in the total light output from a star. Closer home, minute changes in the Earth’s atmosphere and clouds introduce random errors in the observed brightness of the star.

One class of variable stars are eclipsing binaries. When two stars are in orbit around each other, one of the stars sometimes gets eclipsed by the other, as seen from the Earth. This causes the light received at the Earth to drop, and then rise back to the original value as the eclipse ends.

Generally, epoch folding techniques are used to analyse these systems. In this method, a time period is assumed and all observations are wrapped in a circular fashion. For example, if the assumed period is one week, then any observation on day 8 or 15 is the same as day 1. If the assumed time period is correct, the peaks and valleys of the periodic signal will add up and we get a coherent plot, otherwise, we just get a flat profile.

Considering the underlying physics of eclipsing binaries, where the start of an eclipse (a deviation from the mean) has an effect on subsequent observations as the star gets dimmer, it makes sense to model this as some ARMA process with seasonality.

2. Data

The photometry data chosen is of an eclipsing binary observed in the Large Magellanic Cloud (OGLE-III: LMC127.3.10403). There are about 700 observations taken over the course of 4 years.

It has 3 columns: ‘MJD’, ‘i\(\_\)mag’ and ‘err’. MJD stands for Modified Julian Date, a time refernce system similar to UTC used by astronomers. A value ‘MJD = 5260.627’ means the observation was done on day 5260, once 0.627 of the day has passed, i.e. 0.627*24 = 15 hrs. ‘i\(\_\)mag’ is a measure of the observed brightness. ‘err’ is the estimated margin of error in ‘i\(\_\)mag’ reported by the observatory.

3. Visual exploration

The observations are done at non-constant times of the day, as and when telescope time is available. To convert them into equi-spaced observations, aggregate all observations in a 2-day period. Since the period of variability is more than 50 days, it is expected that the brightness will not vary significnatly over two days, and we can thus aggregate the data.

There are 3 spans of continuous observation of \(\sim 100\) data-points that can be used for analysis. Select the observations around MJD = 5900 as this series is the longest. We keep the other two for verification.

The process appears to be mean-stationary and variance stationary. There is a definite seasonal component with period around 30 observation cycles, i.e. 60 days. We first plot the spectral power to determine the period of seasonality.

## [1] "Most dominant frequency, f =  0.0296296296296296  with cycle time  =  33.75 observations."

4. Fitting the ARMA model

4.1 Determining period

There isn’t enough resolution at the lower frequencies to state if the period is 33 or 34 so we fit a SARMA\([1,0]_{t}\) model with \(t\) varying in 30-35.

AIC
Period = 30 samples -448.93
Period = 31 samples -498.09
Period = 32 samples -546.55
Period = 33 samples -645.28
Period = 34 samples -547.10
Period = 35 samples -499.78

4.2 Ascertaining the SARMA model

Now do a sweep of SARMA\([0, 0]\times[p, q]_{33}\) models, with \(p, q \in (0:3)\). We do not take more than 3 lags of period 33, since we have limited data. Therefore, our initial model is of the form \[ y_{n} - \mu = \Phi_{1}(y_{n-33}-\mu) + \Phi_{2}(y_{n-66}-\mu) + \Phi_{3}(y_{n-99}-\mu) + \epsilon_{n} + \Psi_{1}\epsilon_{n-33} + \Psi_{2}\epsilon_{n-66} + \Psi_{3}\epsilon_{n-99}\] where \(y_{i}\)’s are the observations, \(\Phi_{i}, \Psi_{i}\) are the model parameters, \(\mu\) is the intercept term, and \(\epsilon_{i}\)’s follow a Gaussian White-Noise distribution.

Tabulating the AIC values:

SMA0 SMA1 SMA2 SMA3
SAR0 -358.72 -447.93 -493.91 -517.22
SAR1 -645.28 -573.32 -642.68 -571.02
SAR2 -573.01 -572.44 -571.02 -569.02
SAR3 -573.02 -571.02 -569.02 -567.02

The model with best AIC has only SAR1 coefficient.

4.3 Intra seasonal ARMA terms

Now to find the ARMA coefficients. Again, we do not fit beyond 4 lags due to data length constraints, considering the possibility of overfitting.

AIC table:

MA0 MA1 MA2 MA3 MA4
AR0 -645.28 -587.05 -602.92 -610.48 -620.54
AR1 -634.45 -638.46 -637.26 -643.90 -643.48
AR2 -641.95 -641.33 -639.38 -686.75 -640.37
AR3 -643.11 -678.92 -637.37 -684.68 -640.95
AR4 -648.16 -680.32 -675.56 -640.46 -639.02

SARMA\([3, 1]\times[1,0]_{33}\), SARMA\([2, 3]\times[1,0]_{33}\) and SARMA\([4, 1]\times[1,0]_{33}\) models show the best AIC. There is evidence of imperfect maximisation at SARMA\([3, 2]\times[1,0]_{33}\), where the AIC increases by more than 2 on adding a parameter. We prefer smaller models with less coefficients, so the final model to be chosen is SARMA\([3, 1]\times[1,0]_{33}\).
\[ (1-\Phi_{1}B^{33})(1-\phi_{1}B - \phi_{2}B^{2} - \phi_{3}B^{3})(y_{n}-\mu) = (1-\Psi_{1}B^{33})(1-\psi_{1}B )\epsilon_{n}\]

5. Validity of fitted SARMA model

We see that the fitted values are very close to the true data. However, this can simply be due to overfitting as we have used a complex model for a relatively short time-series.

## Series: new_data_3$i_mag 
## ARIMA(3,0,1)(1,0,0)[33] with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ar3      ma1    sar1     mean
##       1.4652  -0.0449  -0.4688  -0.9423  0.0438  17.5997
## s.e.  0.0785   0.1526   0.0786   0.0322  0.1121   0.0021
## 
## sigma^2 estimated as 0.0002984:  log likelihood=346.46
## AIC=-678.92   AICc=-678.01   BIC=-658.8

The SAR1 coefficient of 0.0438 seems rather low for a series that has such strong seasonality

5.1 Root Analysis of ARMA polynomials

## [1] "Absolute values of roots of 3-degree AR polynomial"
## [1] 1.011420 1.011420 2.085124
## [1] "Absolute values of roots of 1-degree MA polynomial"
## [1] 1.061207

The roots of the ARMA polynomial all have moduli greater than 1, indicating a causal, invertible model. But some of the roots are just outside the unit circle, so we need to check whether it is alright to use this model, or reduce some parameters.

5.2 Residual analysis for White-Noise

## [1] "Mean of residuals =  -0.000173687762435632"
## [1] "Median of residuals =  -0.00142491007564629"
## [1] "St. Dev of residuals =  0.0172744771643496"
## [1] "Mean err of raw data =  0.0135572519083969"

## [1]  22 105

There is no significant autocorrelation in the residuals. They have mean and median equal to 0. The Q-Q plot shows that the residuals are almost normally distributed, except for some deviation at the upper tail. Moreover, the predicted value of \(\sigma\) is very close to the average error in the source data. This is consistent with the Gaussian White-Noise assumption.

5.3 Stability of model

Since the three chunks of 100 observations are just different realisations of the same underlying process, we expect that fitting the same SARMA\([3, 1]\times[1,0]_{33}\) model to those data will give similar results.

We cannot directly compare AIC values across the three time series, but we know that, \[\text{Log likelihood, } l = \sum_{i=1}^{n}\log(p_{i})\]

If the true parameter values are the same across time series, then \(p_{i}\) values will be similar. Hence we compare \[\bar l = \sum_{i=1}^{n}\frac{\log(p_{i})}{n}\]

the average log likelihood for each model fit.

dataset number of observations aic avg. log likelihood sar1 intercept sigma MA root moduli AR root moduli
1 131 -678.9235 2.644746 0.0438072 17.59973 0.0172745 1.061207 c(1.01141998560855, 1.01141998560855, 2.08512365744271)
2 117 -573.5340 2.510829 0.5755525 17.60919 0.0189723 1.900688 c(1.13788782032393, 1.95682682354446, 2.53532201973609)
3 105 -466.5663 2.288411 0.3942209 17.58811 0.0245278 2.019870 c(1.22452162629986, 2.37980199623981, 2.37980199623981)

The values of the coefficients differ across realisations. SAR coefficient is larger for datasets 2 and 3, justifying the seasonality assumption. The model still retains causality and invertibility.
But, for dataset-3, the increased white-noise variance and lesser log-likelihood indicate that the model has less predictive power. Therefore, we cannot proceed with this model. If by any means we are able to concatenate the three time series and use the entire history of \(\sim\) 400 observations, perhaps we can get a better estimate.

6. Conclusion

7. References