The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
Getting Started with iAR
Felipe Elorrieta, Cesar Ojeda, Susana Eyheramendy,
Wilfredo Palma
2026-05-16
Overview
The iAR package provides autoregressive models for
time series observed at irregularly spaced time points, a
common situation in astronomy, finance, and environmental science. It
implements three complementary model classes:
iAR |
Scalar \(\phi \in (0,1)\) |
Positive only |
Eyheramendy et al. (2018) |
CiAR |
Vector \((\phi_R, \phi_I)\), \(\|\phi\| < 1\) |
Positive and negative |
Elorrieta et al. (2019) |
BiAR |
Vector \((\phi_R, \phi_I)\), \(\|\phi\| < 1\) |
Bivariate, cross-correlated |
Elorrieta et al. (2021) |
All three share the same workflow: construct →
simulate → estimate → forecast
/ interpolate.
Generating Irregular Observation Times
Real astronomical or environmental data is rarely sampled on a
regular grid. The gentime() function generates realistic
irregular time grids from several distributions. The first argument
x defaults to NULL, so the
utilities object is created internally and you can call
gentime() directly without any setup step.
set.seed(2847)
# x defaults to NULL — no setup needed
times <- gentime(n = 100)@times
cat("Number of observations:", length(times), "\n")
#> Number of observations: 100
cat("Time span: [", round(min(times), 1), ",", round(max(times), 1), "]\n")
#> Time span: [ 162.8 , 3137.7 ]
cat("Mean gap between observations:", round(mean(diff(times)), 2), "\n")
#> Mean gap between observations: 30.05
The default distribution ("expmixture") mixes two
exponential distributions, mimicking data sets where most observations
are closely spaced but occasional long gaps occur. Alternative
distributions include "uniform",
"exponential", and "gamma".
gaps <- diff(times)
hist(log1p(gaps),
main = "Inter-observation gaps (log1p scale)",
xlab = "log(1 + gap)", col = "steelblue", border = "white",
breaks = 20)
Distribution of inter-observation gaps (log
scale) for 100 irregular times drawn from an exponential mixture.
iAR: Univariate Irregular Autoregressive Model
Simulation
Construct an iAR object with the desired coefficient,
then call sim().
set.seed(2847)
times <- gentime(n = 100)@times
# Build the model: family = "norm", phi = 0.9
m_iar <- iAR(family = "norm", times = times, coef = 0.9)
m_iar <- sim(m_iar)
head(m_iar@series)
#> [1] 0.6949442 1.7577536 1.2249041 1.7570633 1.8224373 -1.6668806
plot(m_iar, type = "o",main = "Simulated iAR series (phi = 0.9)",ylab="Values",xlab="Time",pch=20)
Simulated iAR series (phi = 0.9). The irregular
spacing is visible in the unequal distances along the x-axis.
Parameter Estimation
Two estimators are available:
loglik() — direct MLE; supports all
three families ("norm", "t",
"gamma").kalman() — Kalman-filter MLE;
available for all three model classes (iAR,
CiAR, BiAR).
Before estimation, standardise the series (zero mean, unit variance)
and supply a vector of error standard deviations. For simulated data
with no measurement noise, zeros suffice.
y <- m_iar@series
y_std <- y / sd(y) # unit variance
m_iar@series <- y_std
m_iar@series_esd <- rep(0, length(times))
Using loglik()
m_loglik <- loglik(m_iar)
cat("loglik estimate: phi =", round(m_loglik@coef, 4), "\n")
#> loglik estimate: phi = 0.8883
cat("Log-likelihood: ", round(m_loglik@loglik, 4), "\n")
#> Log-likelihood: 111.6524
Using kalman()
m_kalman <- kalman(m_iar)
cat("kalman estimate: phi =", round(m_kalman@coef, 4), "\n")
#> kalman estimate: phi = 0.8876
cat("Kalman log-lik: ", round(m_kalman@kalmanlik, 4), "\n")
#> Kalman log-lik: 0.4045
Both estimators should recover a value close to the true \(\phi = 0.9\).
Standard Errors via the Hessian
Set hessian = TRUE to obtain standard errors and
p-values.
m_h <- iAR(family = "norm", times = times, coef = 0.9, hessian = TRUE)
m_h@series <- y_std
m_h@series_esd <- rep(0, length(times))
m_h <- loglik(m_h)
# Summary contains the coefficient, SE, p-value, information criteria, and residual diagnostics
summary(m_h)
#> iAR model
#>
#> Family:
#> norm
#>
#> Coefficients:
#> Estimate St. Error t value Pr(>|t|)
#> phi 0.89 0.02 38.71 0.00
#>
#> Information criteria:
#> AIC: 225.30
#> BIC: 227.91
#>
#> Residual diagnostics:
#> ACF lag 1: 0.06
#> Ljung-Box test p-value (lag=1): 0.57
#> Ljung-Box test p-value (lag=10): 0.98
Forecasting
After estimation, use forecast() to predict
tAhead time units into the future.
m_fc <- forecast(m_kalman, tAhead = 10)
cat("Forecast (10 time units ahead):", round(m_fc@forecast, 4), "\n")
#> Forecast (10 time units ahead): 0.0761
plot_forecast(m_fc,ylab="Values",xlab="Time",type="o",pch=20)
iAR series with one-step-ahead forecast (blue
dot).
Interpolation (Imputing Missing Values)
interpolation() imputes NA values in the
series using the fitted model.
# Introduce a missing value at position 10
m_interp <- m_kalman
m_interp@series[10] <- NA
m_interp <- interpolation(m_interp)
cat("Interpolated value at position 10:",
round(m_interp@interpolated_values, 4), "\n")
#> Interpolated value at position 10: 1.1321
cat("Original value at position 10 :",
round(y_std[10], 4), "\n")
#> Original value at position 10 : 1.6115
CiAR: Complex Irregular Autoregressive Model
CiAR extends iAR to allow negative
autocorrelation via a complex coefficient \(\phi = \phi_R + i\,\phi_I\). The
stationarity condition is \(|\phi| =
\sqrt{\phi_R^2 + \phi_I^2} < 1\).
Simulation and Estimation
set.seed(2847)
times <- gentime(n = 100)@times
# phi = (0.9, 0): purely real CiAR — same as iAR but more general
m_ciar <- CiAR(times = times, coef = c(0.9, 0))
m_ciar <- sim(m_ciar)
# Standardise
y_c <- m_ciar@series
y_c_std <- y_c / sd(y_c)
m_ciar@series <- y_c_std
m_ciar@series_esd <- rep(0, length(times))
m_ciar <- kalman(m_ciar)
cat("CiAR estimate: phiR =", round(m_ciar@coef[1], 4),
" phiI =", round(m_ciar@coef[2], 4), "\n")
#> CiAR estimate: phiR = 0.9076 phiI = 0
cat("Modulus |phi| =", round(sqrt(sum(m_ciar@coef^2)), 4), "\n")
#> Modulus |phi| = 0.9076
Forecasting
m_ciar <- forecast(m_ciar, tAhead = 10)
cat("CiAR forecast:", round(m_ciar@forecast, 4), "\n")
#> CiAR forecast: 0.8211
BiAR: Bivariate Irregular Autoregressive Model
BiAR models two time series observed at the
same irregular time points with a shared autoregressive
structure and cross-correlation \(\rho\).
Simulation
set.seed(2847)
times <- gentime(n = 100)@times
# coef = c(phiR, phiI), rho = cross-series correlation
m_biar <- BiAR(times = times, coef = c(0.9, 0.3), rho = 0.9)
m_biar <- sim(m_biar)
head(m_biar@series) # matrix: column 1 = series 1, column 2 = series 2
#> [,1] [,2]
#> [1,] 0.9385191 1.4786262
#> [2,] -0.5465484 1.6675406
#> [3,] -0.8640927 0.5525127
#> [4,] -0.9372092 0.3941008
#> [5,] -0.9733738 0.1221019
#> [6,] -0.6999326 -0.9973335
Estimation
n <- length(times)
y_b <- m_biar@series
y_b_std <- y_b / matrix(apply(y_b, 2, sd), nrow = n, ncol = 2, byrow = TRUE)
m_biar@series <- y_b_std
m_biar@series_esd <- matrix(0, n, 2)
m_biar <- kalman(m_biar)
cat("BiAR estimate: phiR =", round(m_biar@coef[1], 4),
" phiI =", round(m_biar@coef[2], 4), "\n")
#> BiAR estimate: phiR = 0.9208 phiI = 0.2989
Forecasting
m_biar <- forecast(m_biar, tAhead = 10)
cat("BiAR forecast (series 1):", round(m_biar@forecast[1], 4), "\n")
#> BiAR forecast (series 1): -0.2926
cat("BiAR forecast (series 2):", round(m_biar@forecast[2], 4), "\n")
#> BiAR forecast (series 2): 0.3773
Real Data Example: AGN Light Curve
The package includes the agn dataset — 237 K-band flux
measurements of the active galactic nucleus MCG-6-30-15, taken between
2006 and 2011.
data(agn)
head(agn)
#> t m merr
#> 1 4087.834 2.878444e-15 6.709391e-17
#> 2 4091.857 2.803222e-15 6.377840e-17
#> 3 4095.849 2.755297e-15 6.256418e-17
#> 4 4099.835 2.759666e-15 6.825679e-17
#> 5 4104.821 2.629979e-15 6.472206e-17
#> 6 4107.828 2.728070e-15 6.211146e-17
plot(agn$t, agn$m, xlab = "Heliocentric JD - 2450000",
ylab = "Flux [10^-15 erg/s/cm^2/A]",
main = "AGN MCG-6-30-15 (K band)",type="o",pch=20)
AGN MCG-6-30-15 light curve. Gaps in the
observation schedule produce the irregular spacing typical of
ground-based astronomical surveys.
Fitting an iAR Model
# Standardise: centre and scale by SD (measurement errors available)
y_agn <- agn$m
y_agn_c <- (y_agn - mean(y_agn)) / sd(y_agn)
esd_agn <- agn$merr / sd(y_agn)
m_agn <- iAR(
family = "norm",
times = agn$t,
series = y_agn_c,
series_esd = esd_agn
)
m_agn <- kalman(m_agn)
#> Warning: Estimated coefficient (0.9973) is near the upper boundary of (0, 1).
#> The model may be near non-stationary; results should be interpreted with
#> caution.
cat("Estimated phi:", round(m_agn@coef, 4), "\n")
#> Estimated phi: 0.9973
cat("Kalman log-lik:", round(m_agn@kalmanlik, 4), "\n")
#> Kalman log-lik: -1.2607
Fitted Values and Forecast
m_agn <- fit(m_agn)
plot_fit(m_agn, xlab = "Heliocentric JD - 2450000",
ylab = "Flux [10^-15 erg/s/cm^2/A]",type="o",pch=20)
AGN light curve (standardised) with fitted iAR
values overlaid.
m_agn <- forecast(m_agn, tAhead = 1)
cat("One-step-ahead forecast:", round(m_agn@forecast, 4), "\n")
#> One-step-ahead forecast: -1.5625
Utility: Pairing Two Irregular Time Series
When two series are observed at different (but overlapping) time
points, pairingits() matches them within a tolerance window
before fitting a BiAR model.
data(cvnovag)
data(cvnovar)
# Pair the G-band and R-band CV Nova light curves
paired <- pairingits(data1 = cvnovag, data2 = cvnovar, tol = 0.01)
cat("Paired observations:", nrow(paired@paired), "\n")
#> Paired observations: 129
head(paired@paired)
#> t m merr
#> [1,] NA NA NA 58278.29 18.32443 0.05345725
#> [2,] NA NA NA 58281.32 18.34270 0.05294770
#> [3,] NA NA NA 58284.33 18.46952 0.06193419
#> [4,] NA NA NA 58287.29 18.41748 0.05424246
#> [5,] 58290.23 18.64738 0.1674221 NA NA NA
#> [6,] 58293.34 18.39692 0.1456636 NA NA NA
Summary
The table below collects the key functions for each step of the
modelling workflow.
| Construct |
iAR(family, times, coef) |
CiAR(times, coef) |
BiAR(times, series, coef, rho) |
| Simulate |
sim(m) |
sim(m) |
sim(m) |
| Estimate (direct MLE) |
loglik(m) |
— |
— |
| Estimate (Kalman) |
kalman(m) |
kalman(m) |
kalman(m) |
| Forecast |
forecast(m, tAhead) |
forecast(m, tAhead) |
forecast(m, tAhead) |
| Interpolate |
interpolation(m) |
interpolation(m) |
interpolation(m) |
| Plot fitted |
plot_fit(m) |
plot_fit(m) |
plot_fit(m) |
| Plot forecast |
plot_forecast(m) |
plot_forecast(m) |
plot_forecast(m) |
| Generate times |
gentime(n = ...) |
— |
— |
| Pair two series |
pairingits(data1, data2) |
— |
— |
References
Eyheramendy, S., Elorrieta, F., & Palma, W. (2018). An irregular
discrete time series model to identify residuals with autocorrelation in
astronomical light curves. Monthly Notices of the Royal Astronomical
Society, 481(4), 4990–5003. doi:10.1093/mnras/sty2487
Elorrieta, F., Eyheramendy, S., & Palma, W. (2019). Discrete
semi-parametric model for the investigation of variable stars.
Astronomy & Astrophysics, 627, A120. doi:10.1051/0004-6361/201935560
Elorrieta, F., Eyheramendy, S., Palma, W., & Ojeda, C. (2021). A
bivariate irregular autoregressive model. Monthly Notices of the
Royal Astronomical Society, 505(1), 1105–1116. doi:10.1093/mnras/stab1216
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.