1 Introduction
Time series analysis is typically performed on data observed at regular intervals, under the assumption of evenly spaced discrete time points (see, for instance, Box and Jenkins (1970); Brockwell and Davis (1991); Hyndman and Athanasopoulos (2021)). When the data are recorded at irregular time intervals, two main strategies are commonly adopted: modeling the series in continuous time, or transforming it to a regularly spaced sequence through interpolation (see Rehfeld et al. (2011)).
Continuous-time models rely on the assumption that consecutive observations are closely spaced. However, this condition is not always met, especially in applications with sparse or unevenly sampled data. The package introduced in this paper addresses this challenge by providing discrete-time models specifically designed for unequally spaced observations.
Although several general-purpose time series packages are available in R, most notably those described in Hyndman et al. (2026) and Hyndman and Khandakar (2008), which provide a comprehensive framework for modeling, forecasting, and evaluating regularly observed time series, there are currently few tools available for modeling irregularly observed time series. For instance, two notable packages that implement Continuous ARMA (CARMA) models are growth and yuima (Brouste et al. 2014). The main distinction is that growth allows incorporating fixed effects into the CARMA specification, whereas this is not supported in yuima. An extension of the CARMA model allowing a Hurst parameter \(H\) to represent long-term memory is implemented in the carfima package.
CARMA models can also be estimated through more general frameworks for stochastic differential equations with irregular time points, such as the sde (McMillen 2008) and ctsmr (Juhl et al. 2016) packages. These packages use maximum likelihood estimation combined with Kalman filtering techniques.
Multivariate irregular time series can be modeled using Rdrw and ctsem, which implement damped random walk processes. Additionally, several other packages in R are designed for specific tasks on irregular data. For instance, the lomb package (Ruf 1999) implements the Lomb–Scargle periodogram, while the RobPer package (Thieler et al. 2016) allows the estimation of its extension, known as the Generalized Lomb–Scargle (GLS) periodogram (Zechmeister and Kürster 2009), both of which are used for period estimation in irregularly observed time series.
Another example is the sour package (Edelson et al. 2017) which allows one to obtain the Discrete Correlation Function (DCF) with which one can obtain an estimator of the autocorrelation function under irregular sampling. Moreover, the package BINCOR (Polanco-Martinez et al. 2019) addresses the problem of estimating the correlation (and cross-correlation) between two series that are unevenly spaced — and possibly sampled at different timepoints — by using a “binned correlation” approach which resamples the irregular series into regular bins (taking into account persistence/autocorrelation) and then computes correlations or cross-correlations.
This paper presents iAR, a new R package that introduces novel discrete-time models to fit irregularly observed time series. Unlike the packages mentioned above, which primarily rely on continuous-time formulations, iAR adopts a discrete-time perspective, offering greater flexibility in the modeling and estimation process.
The models implemented in iAR include the irregular autoregressive (iAR) process, its complex-valued extension (CiAR), and a bivariate version (BiAR). Additionally, the package provides tools to handle non-Gaussian data through extensions of the iAR model for \(t\) and Gamma-distributed observations. These models allow for the estimation of negative autocorrelation, which cannot be captured by CARMA-based approaches.
iAR offers a range of functionalities for irregularly spaced data, including: (i) random generation of observation times; (ii) simulation of the irregular autoregressive processes; (iii) estimation of autocorrelation; (iv) model fitting via maximum likelihood; and (v) interpolation of missing values and forecasting. Furthermore, all models implemented in iAR are also compatible with regularly spaced stationary series, such as standard AR processes.
The package is built using the S7 object-oriented system (Vaughan et al. 2024), which provides a robust and modular framework for class and method definitions. Compared to earlier systems like S3 and S4, S7 enforces a clearer separation between class specification and method implementation, supports formal validation, and ensures consistent method dispatch. This structure improves code clarity, extensibility, and integration with other object-oriented tools in R. Internally, the package relies on the zoo package (Zeileis and Grothendieck 2005) to manage indexed irregular time series efficiently, supporting multiple time representations, including numeric time scales,
calendar-based indices (e.g., Date and POSIXct), and user-defined irregular sampling schemes.
The remainder of this paper is organized as follows: Section 2 introduces the modeling framework for irregularly observed autoregressive processes and discusses their applicability. Section 3 presents the procedures to obtain the maximum likelihood estimators for each of the models available in the iAR package. Section 4 describes the structure and methods implemented in the package through examples using simulated data. Section 5 reports the computational times of the estimation procedures and compares them with similar functions available in other packages. Section 6 presents applications of the iAR package to astronomical data. Finally, Section 7 concludes the paper with a discussion.
2 Irregularly observed autoregressive models
The five models available in the iAR package are presented in this section and are depicted in Figure 1. The choice of each model depends on the type of time series to be fitted. For example, if the data available are non-Gaussian distributed or not positively autocorrelated the appropriate model to fit the data may vary. In the following, each of the five models is briefly described.
Figure 1: Flowchart of Irregularly Observed Autoregressive Processes.
2.1 iAR Model
The irregular Autoregressive (iAR) model was introduced by Eyheramendy et al. (2018). We present here a brief description of the model. Let \(\{t_j\}\) be the irregular observational times for \(j=1,\ldots,n\). It is assumed that in general \(t_j-t_{j-1}\) are not equal for \(j=2,\ldots,n\). The iAR model is defined by,
\[\begin{equation} y_{t_j}=\phi^{t_j-t_{j-1}} \, y_{t_{j-1}} + \sigma \, \sqrt{1-\phi^{2(t_j-t_{j-1})}} \, \varepsilon_{t_j}, \tag{1} \end{equation}\]
where \(\varepsilon_{t_j}\) is a white noise sequence with zero mean and unit variance, \(\sigma\) is the standard deviation of the sequence of observations \(y_{t_j}\). Note that \(\phi\) is the parameter that describes the autocorrelation of the process. The iAR model can be represented in a state-space system as follows,
\[\begin{equation} X_{t_j} = \phi^{t_j-t_{j-1}} X_{t_{j-1}} + \epsilon_{t_j} \tag{2} \end{equation}\] \[\begin{equation} y_{t_j} = X_{t_j} + \omega_{y_{t_j}} \tag{3} \end{equation}\]
Note that, in this representation the transition matrix of dimension \(1\times1\) is defined as \(F_{t_j} = \phi^{t_j-t_{j-1}}\). On the other hand, the observation matrix is defined by \(G_{t_j} = 1\). The measurement error of \(y_{t_j}\) is defined by \(\omega_{y_{t_j}}\) with known variance \(\delta^2_{t_j}\) and the error of the transition equation is \(\epsilon_{t_j}\) with variance \(Q_{t_j} =\sigma^2(1-\phi^{2(t_j-t_{j-1})}\)).
Assuming Gaussian errors \(\varepsilon_{t_j}\), i.e. \(\varepsilon_{t_j}\sim N(0,1)\), the iAR model is equivalent to the CARMA(1,0) or CAR(1) model. However, the iAR model can be extended to allow for non-Gaussian distributed data. The extensions implemented in the iAR package assume that the conditional moments of the iAR process follow either a Gamma or Student-t distribution. These two models are described next.
2.2 iAR-Gamma Model
The conditional mean and variance of the iAR model are:
\[\begin{eqnarray*} \mathbb{E}(y_{t_j}|y_{t_{j-1}})&=& \mu + \phi^{t_j-t_{j-1}} \, y_{t_{j-1}}\\ \mathbb{V}(y_{t_j}|y_{t_{j-1}})&=& \sigma^2 \left(1-\phi^{2(t_j-t_{j-1})}\right)\\ \end{eqnarray*}\]
respectively. If \(y_{t_j}|y_{t_{j-1}} \sim\) Gamma\((\alpha_{t_j},\beta_{t_j})\), where \(\alpha_{t_j}\) and \(\beta_{t_j}\) are the shape and scale parameters respectively, then \(y_{t_j}\) will be an iAR-Gamma process if \(\alpha_{t_j}\) and \(\beta_{t_j}\) are such that the following two equations are satisfied,
\[\begin{equation} \begin{split} \alpha_{t_j}\beta_{t_j} &= \mu + \phi^{t_j-t_{j-1}} \, y_{t_{j-1}}\\ \alpha_{t_j}\beta_{t_j}^2 &= \sigma^2 \left(1-\phi^{2(t_j-t_{j-1})}\right)\\ \end{split} \tag{4} \end{equation}\]
2.3 iAR-T Model
A similar approach can be used to define an iAR-T model. Suppose that \(y_{t_j}|y_{t_{j-1}} \sim t_{\nu}(\lambda_{t_j},\tau_{t_j}^2)\) is a non-standardized Student-T distribution with mean \(\lambda_{t_j}\), variance \(\tau_{t_j}^2\) and degrees of freedom \(\nu\). Then \(y_{t_j}\) will be an iAR-T process if \(\lambda_{t_j}\) and \(\tau_{t_j}^2\) are defined by,
\[\begin{equation} \begin{split} \lambda_{t_j} &= \mu + \phi^{t_j-t_{j-1}} \, y_{t_{j-1}}\\ \tau_{t_j}^2 &= \frac{\nu-2}{\nu}\sigma^2 \left(1-\phi^{2(t_j-t_{j-1})}\right).\\ \end{split} \tag{5} \end{equation}\]
2.4 CiAR Model
For irregularly observed time series, the models described above are well defined for positive values of the autocorrelation parameter \(\phi\). The reason for this is that the time gap \(\Delta_j = t_j-t_{j-1}\) can take any positive real value. Consequently, the \(\phi\) parameter must be greater than zero since the result of a negative \(\phi\) to the power of a real number, in general, is a complex number.
In order to address this issue, equation (1) can be extended to a complex process which allows for negative values of the autocorrelation of the process. This extension of the iAR model is called the Complex irregular Autoregressive (CiAR) model (Elorrieta et al. 2019). The CiAR process is characterized by a complex autocorrelation parameter defined by \(\phi= \phi^R + i \phi^I\). The CiAR model can be represented by the following state space system,
\[\begin{equation} \left(\begin{array}{c} y_{t_j}^{*} \\ z_{t_j}^{*} \end{array} \right)= \left(\begin{array}{cc} \alpha_{t_j} & -\beta_{t_j} \\ \beta_{t_j} & \alpha_{t_j} \end{array} \right)\left(\begin{array}{c} y_{t_{j-1}}^{*} \\ z_{t_{j-1}}^{*} \end{array} \right) + \sigma \sqrt{1-\left|\phi^{t_j-t_{j-1}}\right|^2} \left(\begin{array}{c} \varepsilon_{y_{t_j}} \\ \varepsilon_{z_{t_j}} \end{array} \right) \tag{6} \end{equation}\] \[\begin{equation} y_{t_j} = \left(\begin{array}{cc} 1 & 0 \end{array} \right) \left(\begin{array}{c} y_{t_j}^{*} \\ z_{t_j}^{*} \end{array} \right) + \omega_{y_{t_j}} \tag{7} \end{equation}\]
In this case, the transition equation is two-dimensional with a transition matrix \(F_{t_j} = \left(\begin{array}{cc} \alpha_{t_j} & -\beta_{t_j} \\ \beta_{t_j} & \alpha_{t_j} \end{array} \right)\). In addition, the observation matrix is defined by \(G_{t_j}=\left(\begin{array}{cc} 1 & 0 \end{array} \right)\). \(\alpha_{t_j}\) and \(\beta_{t_j}\) are the reparameterized coefficients of the model, given by \(\alpha_{t_j} = |\phi|^{\delta_{j}} \cos(\delta_{j} \psi)\) and \(\beta_{t_j} = |\phi|^{\delta_{j}} \sin(\delta_{j} \psi)\) respectively, where \(\delta_{j} = t_j-t_{j-1}\) and \(\psi = \arccos\left( \frac{\phi^{R}}{|\phi|}\right)\). The variance of the errors of the transition and observations equation are defined by \(Q_{t_j} = \sigma^2(1-\left|\phi^{t_j-t_{j-1}}\right|^2)\) and \(\delta^2_{t_j}\) respectively. As in the iAR model, \(\sigma^2\) denotes the marginal variance of the process.
2.5 BiAR Model
We propose a bivariate extension of the CiAR model which we called the bivariate irregular autoregressive (BiAR) (Elorrieta et al. 2021) model. Unlike the CiAR model, the BiAR model is defined by a two-dimensional observation equation. The state space representation of the BiAR model is defined as,
\[\begin{equation} \left(\begin{array}{c} y_{t_j}^{*} \\ z_{t_j}^{*} \end{array} \right)= \left(\begin{array}{cc} \alpha_{t_j} & -\beta_{t_j} \\ \beta_{t_j} & \alpha_{t_j} \end{array} \right)\left(\begin{array}{c} y_{t_{j-1}}^{*} \\ z_{t_{j-1}}^{*} \end{array} \right) + \left(\begin{array}{c} \varepsilon_{y_{t_j}} \\ \varepsilon_{z_{t_j}} \end{array} \right) \tag{8} \end{equation}\] \[\begin{equation} \left(\begin{array}{c} y_{t_j} \\ z_{t_j} \end{array} \right) = \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right) \left(\begin{array}{c} y_{t_j}^{*} \\ z_{t_j}^{*} \end{array} \right) + \left(\begin{array}{c} \omega_{y_{t_j}} \\ \omega_{z_{t_j}} \end{array} \right) \tag{9} \end{equation}\]
where \(\varepsilon_{t_{j}} = \left(\begin{array}{c} \varepsilon_{y_{t_j}} \\ \varepsilon_{z_{t_j}} \end{array} \right)\) can be decomposed as \(\varepsilon_{t_{j}} = L_{t_{j}} \xi_{t_j}\), where \(L_{t_j}\) is a lower triangular matrix obtained by the Cholesky decomposition of the covariance matrix of \(\varepsilon_{t_{j}}\) and \(\xi_{t_j}\) is a bivariate normal random variable with zero mean and variance \(\Sigma_{\xi} = \left(\begin{array}{cc}1 & \rho_{\xi} \\ \rho_{\xi} & 1 \end{array} \right)\). As in the CiAR model, \(\alpha_{t_j}\) and \(\beta_{t_j}\) are a reparameterization of the parameters \(\phi^R\) and \(\phi^I\). Note that, the BiAR model is defined by three parameters \(\phi^R\),\(\phi^I\) and \(\rho_{\xi}\) which can be interpreted as the autocorrelation parameter \(\phi^R\), the cross correlation parameter \(\phi^I\) and the contemporary correlation parameter \(\rho_{\xi}\).
2.6 Considerations
All irregularly observed time series models implemented in this package depend on the observational gaps through both the autoregressive coefficient and the innovation variance. This may lead to some edge cases. For instance, when two consecutive observation times are identical (\(\Delta_j = 0\)), the autoregressive coefficient satisfies \(\phi^{\Delta_j}=1\) and the innovation variance becomes zero. In practice, this may lead to numerical instability during estimation procedures. On the other hand, for very large time gaps, the autoregressive coefficient tends to zero, implying that consecutive observations become approximately independent.
3 Estimation algorithms
3.1 Likelihood of the processes
Assume that the \(y_{t_j}\) sequence follows an iAR model, if we further assume a Gaussian distribution, the negative log-likelihood of the \(y_{t_j}\) sequence is given by,
\[\begin{equation} \ell(\theta)=\frac{n}{2}\log (2\pi)+\frac{1}{2}\sum_{j=1}^n \log \nu_{t_j} + \frac{1}{2}\sum_{j=1}^n \frac{e_{t_j}^2}{\nu_{t_j}}, \tag{10} \end{equation}\]
where \(\theta = (\phi, \sigma^2)\) is the vector of the parameters of the model. The initial values are defined by \(e_{t_1}=y_{t_1}\), \(\nu_{t_1}=\sigma^2 + \delta_{t_1}^2\) and \(\hat{y}_{t_1}=0\), while \(\hat{y}_{t_j} = \hat{\phi}^{t_j-t_{j-1}}\, y_{t_{j-1}}\), \(e_{t_j}=y_{t_j}-\hat{y}_{t_j}\) and \(\nu_{t_j}=\sigma^2 (1-\phi^{2(t_j-t_{j-1})}) + \delta_{t_j}^2\) for \(j=2,\ldots,n\).
For the iAR-Gamma model, the log-likelihood of the conditional distribution of \(y_{t_j}|y_{t_{j-1}}\) can be written as,
\[\begin{equation} \ell(\theta) = - \left(\alpha_{t_j}\right) \log \beta_{t_j} - \log \Gamma \left(\alpha_{t_j}\right) - \frac{1}{\beta_{t_j}}y_{t_j} + \left(\alpha_{t_j}-1\right) \log y_{t_j} \tag{11} \end{equation}\]
where \(\theta = (\alpha_{t_j}, \beta_{t_j})\) is the vector of the parameters of the model. Note from equation (4) that \(\alpha_{t_j}\) and \(\beta_{t_j}\) depend on the parameters \(\phi\), \(\sigma^2\) and \(\mu\). The initial value is defined as having a Gamma distribution with shape and scale parameters equal to one, i.e. \(y_{t_1} \sim {\rm Gamma}(1,1)\).
For the iAR-T model, the log-likelihood of the conditional distribution of \(y_{t_j}|y_{t_{j-1}}\) can be written as,
\[\begin{equation} \ell(\theta) = - \frac{1}{2} \log \tau_{t_j}^2 -\frac{\nu+1}{2} \log \left( 1 + \frac{1}{\nu} \left(\frac{y_{t_j}-\lambda_{t_j}}{\tau_{t_j}} \right)^2 \right) \tag{12} \end{equation}\]
where \(\theta = (\nu, \tau_{t_j}, \lambda_{t_j})\) is the vector of the parameters of the model. Note from equation (5) that \(\tau_{t_j}\) and \(\lambda_{t_j}\) depend on the parameters \(\phi\), \(\sigma^2\). The initial value is defined by \(y_{t_1} \sim N(0,1)\).
3.2 Kalman recursions
The state-space representations of the models described in Section 2 allow one to obtain the maximum likelihood estimators of the parameters of the models using the Kalman recursions. Given the transition matrix \(F_{t_j}\), the observation matrix \(G_{t_j}\), the transition error variance \(Q_{t_j}\) and the observation error variance \(R_{t_j}\) the Kalman recursions can be obtained as,
\[\begin{equation} \begin{split} \Lambda_{t_j} &= G \Omega_{t_j} G' + R_{t_j} \\ \Theta_{t_j} &= F_{t_j}\Omega_{t_j} G' \\ \Omega_{t_{j+1}} &= F_{t_j} \Omega_{t_j} F_{t_j}' + Q_{t_j} - \Theta_{t_j} \Lambda_{t_j}^{-1} \Theta_{t_j}' \\ \nu_{t_j} &= X_{t_j} - G \hat{X}_{t_j}\\ \hat{X}_{t_{j+1}} &= F_{t_j} \hat{X}_{t_j} + \Theta_{t_j} \Lambda_{t_j}^{-1} \nu_{t_j} \end{split} \tag{13} \end{equation}\]
Using the Kalman recursions parameters, we can obtain the maximum likelihood estimator from the log-likelihood function defined as,
\[\begin{equation} \ell(\phi) \propto - \frac{1}{2} \mathop{\sum}\limits_{j=1}^n \left( \log (|\Lambda_{t_j}|) + \nu_{t_j}' \Lambda_{t_j}^{-1} \nu_{t_j} \right) \end{equation}\]
4 Structure and usage of the iAR package
4.1 Object-oriented system
The iAR package is structured using the S7 object-oriented system, which provides a formal and modular framework for defining classes and methods. This design enables a coherent and consistent internal organization across the different autoregressive models implemented in the package: iAR, CiAR, and BiAR.
Each model class provides a standardized set of methods for tasks such as maximum likelihood estimation, computation of fitted values, forecasting, and interpolation. The implementation of S7 ensures that these methods are defined uniformly across models, while still allowing for model-specific behavior when needed. Also, extensive validations have been implemented within the package functions to ensure robust behavior under a variety of input scenarios and model configurations.
For handling irregular time indices, the package leverages the zoo package, which offers flexible and efficient tools for managing time-based data in R. This integration allows for seamless manipulation and alignment of time series with irregular observational times, a fundamental feature for the models supported by the package.
4.2 The R interface
Our package can be installed and loaded in R by the usual commands:
install.packages("iAR")
library("iAR")The classes and methods implemented in the iAR package are summarized in Table1. Each class—iAR, CiAR, and BiAR—includes a core set of methods for simulating time series (sim), estimating parameters via maximum likelihood (kalman), computing fitted values (fit), summarizing fitted models (summary), interpolating missing observations (interpolation), forecasting future values (forecast), and visualizing results through various plotting functions (plot, plot_fit, plot_forecast). The iAR class additionally includes loglik for computing the log-likelihood directly.
| Class | Method | Description |
|---|---|---|
| iAR | sim | Sampling from an iAR model |
| loglik | MLE estimation of iAR model | |
| kalman | MLE estimation of iAR model | |
| fit | Fitted values of iAR model | |
| summary | Summary of the iAR model | |
| interpolation | Interpolation from the iAR model | |
| forecast | Forecast from iAR model | |
| plot | Plot of the series | |
| plot_fit | Plot of the fitted values | |
| plot_forecast | Plot of the forecasted values | |
| CiAR | sim | Sampling from a CiAR model |
| kalman | MLE estimation of CiAR model | |
| fit | Fitted values of CiAR model | |
| summary | Summary of the CiAR model | |
| interpolation | Interpolation from the CiAR model | |
| forecast | Forecast from CiAR model | |
| plot | Plot of the series | |
| plot_fit | Plot of the fitted values | |
| plot_forecast | Plot of the forecasted values | |
| BiAR | sim | Sampling from a BiAR model |
| kalman | MLE estimation of BiAR model | |
| fit | Fitted values of BiAR model | |
| summary | Summary of the BiAR model | |
| interpolation | Interpolation from the BiAR model | |
| forecast | Forecast from BiAR model | |
| plot | Plot of the series | |
| plot_fit | Plot of the fitted values | |
| plot_forecast | Plot of the forecasted values | |
| utilities | gentime | Generating time points based on a statistical distribution |
| pairingits | Pairing the observational times of two time series | |
| phase | Computing phased time series | |
| harmonicfit | Harmonic fit to time series |
Most methods are consistently available across the three main classes, with the only exception being the loglik method. This method is implemented exclusively for the iAR class, while parameter estimation for the CiAR and BiAR classes relies entirely on the Kalman filtering framework. Furthermore, since the same methods are used for the Gaussian and non-Gaussian versions of the iAR model (e.g., iAR-Gamma and iAR-T), the underlying distributional assumptions are controlled via arguments passed to these functions. This unified interface simplifies the application of various models to irregularly observed time series.
Additionally, the utilities class provides support functions that enhance the modeling workflow. The gentime method generates irregular time points based on a statistical distribution, which is useful when the time slot is not provided by the user. The pairingits method allows pairing observational times from two irregularly observed time series, facilitating multivariate modeling with the BiAR class. The phase method computes phased time series, which is particularly useful for visualizing periodic behaviors in irregularly observed data. Finally, the harmonicfit method performs harmonic fitting on time series, providing a flexible tool for removing dominant periodic components before autoregressive modeling.
To improve computational efficiency, several functions in the methods were implemented in C++ using the package Rcpp (Eddelbuettel and Francois 2011). In particular, we use the C++ linear algebra library Armadillo, accessed via the RcppArmadillo package (Eddelbuettel and Sanderson 2014).
4.3 Examples
The sim method allows the user to generate realizations from any of the five models available in this package. To simulate a standard iAR process, we first create an object of class iAR as follows:
model_norm <- iAR(family = "norm", coef = 0.9) Here, the model is defined by specifying the distribution family and the autoregressive coefficient. To simulate an iAR process, the argument coef, which corresponds to the model autoregressive parameter, must take a value between 0 and 1.
The time points can be specified by the user when the iAR object is defined. However, if the times argument is left empty, the sim method will automatically generate a time array of length n, using a random distribution defined by the gentime method from the utilities class. By default, an array of size n=100 is generated from a mixture of two exponential distributions with rate parameters lambda1 = 130 and lambda2 = 6.5, and weights p1 = 0.15 and p2 = 0.85, respectively.
Once the model is defined, the process can be simulated using the sim method as follows:
model_norm <- sim(model_norm) To visualize the simulated process, the plot method from the iAR package can be used:
Figure 2: Simulated Irregularly Observed Autoregressive process.
The simulated iAR process is shown in Figure 2. In this figure, the irregular time gaps and the stationary behavior of the process can be observed. To estimate the parameter phi, the kalman method is available in our package. This method uses Kalman recursions to obtain the maximum likelihood estimator based on the state-space representation of the iAR model. Applying this method to the previously simulated process yields the following result:
Another option for estimating the phi parameter of the iAR model is to use the loglik method. This method computes the maximum likelihood estimator based on the full likelihood function. The loglik method can be applied as follows:
To assess whether the parameter \(\phi\) is significantly different from zero, the summary method can be used. This method returns the p-value of the t-test associated with this parameter. To use it properly, the hessian slot must be set to TRUE in the iAR object beforehand. If the p-value returned by this function is less than a given significance level, we can conclude that the parameter is significantly greater than zero. Furthermore, the summary method also provides additional information about the fitted model, including information criteria and residual diagnostics. The usage of this method is as follows:
iAR model
Family:
norm
Coefficients:
Estimate St. Error t value Pr(>|t|)
phi 0.85 0.03 29.51 0.00
Information criteria:
AIC: 237.51
BIC: 240.11
Residual diagnostics:
ACF lag 1: -0.04
Ljung-Box test p-value (lag=1): 0.68
Ljung-Box test p-value (lag=10): 0.71 As shown in this example, the obtained p-value is very small. Consequently, the parameter is significantly greater than zero. Furthermore, the diagnostics of the fitted model suggest no remaining autocorrelation in the residuals, supporting the whiteness of the errors.
For the remaining models available in this package, the methods are structured similarly to the previous example (as shown in Table 1). Thus, the procedure for generating a non-Gaussian iAR model begins by creating an object of the iAR object with either a gamma or t distribution. Below, we illustrate the simulation of an iAR-Gamma process.
where the additional parameters mean and variance, representing the mean and variance of the simulated process, respectively, must be defined. Figure 3 shows a simulated iAR-Gamma process with the distribution of the simulated data.
par(mfrow=c(1,2))
plot(model_gamma,main="(a)",ylab="Values",xlab="Time",type = "o",pch=20)
hist(model_gamma@series,breaks=25,main="(b)",xlab="Values")
Figure 3: Simulated iAR-Gamma process. Figure a) shows the simulated irregularly observed time series. Figure b) shows the histogram of the values of the simulated process.
Note in Figure 3-b that the simulated values have a left-skewed distribution typical of the gamma family. To estimate the parameters of the iAR-Gamma process, the loglik method can be used as follows:
iAR model
Family:
gamma
Coefficients:
Estimate St. Error t value Pr(>|t|)
phi 0.90 0.01 83.80 0.00
mean 1.01 0.04 22.75 0.00
variance 0.83 0.14 5.88 0.00
Information criteria:
AIC: 179.25
BIC: 187.07 In this example, the estimated values for the three parameters are close to their true values, indicating that the estimation procedure performs well. Additionally, the near-zero p-values confirm the statistical significance of the parameter estimates. Note that we use diagnostic = FALSE in the summary method, which suppresses the residual diagnostics to simplify the presentation of the output.
To generate a heavy-tailed iAR process, the method sim can be applied to an iAR object of the family t. The family t requires the additional argument df which specifies the degrees of freedom of the conditional Student’s t distribution. The procedure to generate an iAR-t process is as follows:
par(mfrow=c(1,2))
plot(model_t,main="(a)",ylab="Values",xlab="Time",type = "o",pch=20)
hist(model_t@series,breaks=25,main="(b)",xlab="Values")
Figure 4: Simulated iAR-T process. Figure a) shows the simulated irregularly observed time series. Figure b) shows the histogram of the values of the simulated process.
Note in Figure 4 that the simulated process exhibits a heavy tail distribution. In order to estimate the parameters of this process via maximum likelihood, the method loglik must be used as follows:
iAR model
Family:
t
d.f.:
3
Coefficients:
Estimate St. Error t value Pr(>|t|)
phi 0.90 0.02 42.62 0.00
sigma 0.96 0.23 4.25 0.00
Information criteria:
AIC: 175.25
BIC: 180.46 This method allows estimation of the parameters phi and sigma from the iAR-T model. The standard errors and p-values indicate that the estimates are statistically significant.
To generate a CiAR process, a CiAR class object must first be created. The argument coef, which corresponds to a vector containing the coefficients of the CiAR model must be defined. Finally, the sim method is used to simulate the CiAR process. Note that in the first examples, the observational times are not provided, and are therefore generated by default. To illustrate how observational times can be specified by the user, the gentime function is manually used to generate them. In this case, a mixture of two exponential distributions is considered with means lambda1 = 15 and lambda2 = 2.5, and weights p1 = 0.15 and p2 = 0.85, respectively, as follows:
To estimate the parameters of the CiAR model, the kalman method must be used. If the argument hessian = TRUE is specified, the significance of the estimated parameters can also be obtained, as shown below:
CiAR model
Coefficients:
Estimate St. Error t value Pr(>|t|)
phiR 0.91 0.10 8.86 0.00
phiI -0.04 0.48 -0.08 0.93
Information criteria:
AIC: 298.27
BIC: 304.86 As shown, this method provides both the parameter estimates and the reduced negative log-likelihood evaluated at the coefficients phiR and phiI. The estimates obtained are accurate for both parameters.
Additionally, to obtain the fitted values from the CiAR model, the fit method must be applied. This method returns the fitted values of a time series observed at irregular time points, given the parameters phiR and phiI, as follows:
Figure 5: Simulated Complex Irregularly Observed Autoregressive process. The blue line corresponds to the fitted values estimated with the Kalman recursions.
Figure 5 shows the fitted values estimated with the method fit over the simulated CiAR process. The method forecast applied on a CiAR object allows forecasting from the CiAR model given the coefficients and the tAhead which corresponds to the time gap in which the forecast is made. In the following example, a forecast 1.5 days in the future is performed:
forCIAR <- forecast(model_CiAR, tAhead=1.5)
forCIAR@forecast[1] -0.6434017It is important to notice that the forecast predictor of the CiAR process (as well as for any autoregressive process) converges to the mean for large prediction horizons.
The sim method applied on a CiAR object also allows one to generate an irregularly observed autoregressive process with negative autocorrelation as follows:
par(mfrow=c(1,2))
model_CiAR <- CiAR(times = times@times,coef = c(0.9, 0))
model_CiAR <- sim(model_CiAR)
plot(model_CiAR,ylab="Values",xlab="Time",main="(a)",type = "o",pch=20)
model_CiAR <- CiAR(times = times@times,coef = c(-0.9, 0))
model_CiAR <- sim(model_CiAR)
plot(model_CiAR,ylab="Values",xlab="Time",main="(b)",type = "o",pch=20)
Figure 6: Simulated Complex Irregularly Observed Autoregressive processes with positive autocorrelation of order one (Figure (a)) and negative autocorrelation of order one (Figure (b)).
Figure 6 shows the CiAR process simulated with both positive and negative autocorrelation. Note that the time series with negative autocorrelation shows more frequent oscillations around zero, while the positively correlated process exhibits a more persistent behavior.
The last model implemented in this package is the BiAR process. The BiAR class provides functionality for working with this model. To simulate a BiAR process, the argument coef must be defined as a vector of length two, corresponding to the model coefficients. As we mentioned for the previous examples, the time points can be provided by the user or automatically generated by a random distribution. Additionally, a parameter rho can be specified to indicate the correlation between the two components of the bivariate process. The use of the sim method in this class is as follows:
par(mfrow=c(2,1))
plot(model_BiAR@times,model_BiAR@series[,1],main="(a)",xlab="Time",ylab="Values",
type = "o",pch=20)
plot(model_BiAR@times,model_BiAR@series[,2],main="(b)",xlab="Time",ylab="Values",
type = "o",pch=20)
Figure 7: Simulated Bivariate Irregularly Observed Autoregressive process.
Figure 7 shows the simulated BiAR process. As expected, the two time series display similar behavior due to the value of the contemporary correlation parameter \(\rho=0.9\). To estimate the parameters of the BiAR model and assess their significance, the method kalman should be used with the argument hessian = TRUE, as shown below:
model_BiAR@hessian=TRUE
model_BiAR <- kalman(model_BiAR, niter = 20)
summary(model_BiAR,diagnostic=FALSE)BiAR model
Coefficients:
Estimate St. Error t value Pr(>|t|)
phiR 0.90 0.01 139.46 0.00
phiI 0.31 0.01 25.49 0.00
Information criteria:
AIC: -3902.02
BIC: -3896.81 In addition, the method fit allows estimation of the contemporary correlation parameter rho. The usage of this method is the following:
model_BiAR <- fit(model_BiAR)
model_BiAR@rho[1] 0.9232303These results show accurate estimation of the three parameters of the BiAR model.
5 CPU Performance
In this section, we evaluate the computational performance of the functions implemented in the iAR package for R. First, the computational times are compared with respect to other functions available in R for fitting irregularly spaced time series models. The functions included in the comparison are the following: carma function from the growth package, drw function from the Rdrw package and the loglik method from the iAR class of our iAR package.
The microbenchmark function from the package of the same name is used for the comparison of the computational times of the functions mentioned above. For this experiment, an irregularly observed time series of length 200 was generated. Each function was executed 100 times on the generated time series. The distribution of execution times (in microseconds) is presented in Figure 8. Note that the loglik method from our package estimates the model parameters in significantly less time than the other functions analyzed.
Figure 8: Distribution of the computation times (in microseconds) of the functions carma (growth package), drw (Rdrw package) and the loglik method (iAR package).
In a second experiment, we assessed how the computational time of the estimation methods implemented in our package scales with the sample size of simulated time series. Figure 9 shows the mean computation times of the and methods from the class, and the method from the and classes, for sample sizes ranging from 50 to 1000 observations, after 100 repetitions. The results show an approximately linear growth of computation time with sample size. As expected, the estimation methods for the model are considerably faster than those for the and models, reflecting their lower computational complexity.
Figure 9: Mean computation time (in seconds) of the estimation methods in the iAR package.
In addition, Table 2 reports the mean CPU times (in seconds). For example, the iAR model increases from 4^{-4} seconds (n = 50) to 0.0032 seconds (n = 1000), while the CiAR and BiAR models increase from 0.0076 to 0.1612 and from 0.0107 to 0.2163 seconds, respectively. These results are consistent with the approximately linear growth observed in the previous figure. In particular, when n increases by a factor of 20, the computation time increases by a factor of about 20 for CiAR and BiAR, and by a factor of about 5 for iAR.
| n | iAR | CiAR | BiAR |
|---|---|---|---|
| 50 | 0.0004 | 0.0076 | 0.0107 |
| 100 | 0.0006 | 0.0144 | 0.0244 |
| 250 | 0.0011 | 0.0373 | 0.0594 |
| 500 | 0.0018 | 0.0728 | 0.1154 |
| 1000 | 0.0032 | 0.1612 | 0.2163 |
6 Applications of the iAR package on real datasets
6.1 Identification of Multiperiodic time series
In astronomy, there are surveys where finding irregular time series with periodical behavior is not rare. Furthermore, some astronomical objects can have more than one significant period. Generally, this happens in specific types of variable stars, which are called multi-periodic variable stars. In the iAR package an example of this type of irregular time series can be found. This dataset is called dmcep and corresponds to the time series of the magnitude of the brightness of a double mode cepheid (DMCEP) variable star observed by the OGLE survey (Udalski et al. 2015). This time series has been irregularly observed and each observation has associated an observational error. The dataset contains 191 observations and three attributes: the observational times (t), the brightness magnitudes (m), and the known variance of the measurement errors (merr). The dominant frequency of this time series is 0.741, while the second frequency is 0.543. Figure 10 (a) shows the raw time series of this DMCEP variable star. In addition, we folded this time series by its first period, using the method phase of this package. The folded time series is presented in the Figure 10 (b) and reflects the periodical behavior of the brightness of this astronomical object.
Due to the periodic behavior of this irregular time series, a harmonic model must be fitted to it. For this purpose the harmonicfit method of the iAR package can be used. This function requires specifying the number of harmonic components (nham) to include in the model. In this application, the default value nham=4 was used. Later, we fitted an iAR model to the residuals of the harmonic model to verify if the whole time dependence structure of this time series was explained by this model. The parameter estimated by the iAR model was \(\hat{\phi} = 0.02\) , with an associated p-value of 0.42, suggesting no significant autocorrelation after the harmonic fit.
data(dmcep)
f1=0.741
o1<-harmonicfit(data=dmcep,f1=f1)
y=o1@residuals
y=y/sd(y)
merr=dmcep$merr/sd(y)
x <- iAR(family = "norm", times = dmcep[,1], series = y, series_esd = merr, hessian=TRUE)
x = loglik(x)
summary(x)iAR model
Family:
norm
Coefficients:
Estimate St. Error t value Pr(>|t|)
phi 0.02 0.02 0.81 0.42
Information criteria:
AIC: 528.77
BIC: 532.02
Residual diagnostics:
ACF lag 1: -0.13
Ljung-Box test p-value (lag=1): 0.06
Ljung-Box test p-value (lag=10): 0.03 However, the residual diagnostics indicate some remaining dependence structure. This result is further confirmed by fitting the CiAR model to the residuals of the harmonic model, which yields an estimated autocorrelation parameter of \(-0.63\). This result indicates the presence of negative autocorrelation in the residuals, which explains why it is captured by the CiAR model.
Figure 10 (d) shows that the residuals of the harmonic fit folded using the second period of this variable star exhibit a clear serial correlation. Later, this serial correlation can be modeled using a new harmonic model that considers the second period.
Figure 10: Figure a) shows the raw DMCEP time series and Figure b) shows this time series folded by the first period of the variable star. Figure c) shows the residuals of the harmonic fit and Figure d) shows this residuals folded by the second period of the variable star.
6.2 Pairing two astronomical time series
Among the models available in the iAR package, the BiAR model allows two time series to be modeled jointly. However, this model assumes that two time series are measured at the same times. This assumption is not necessarily satisfied. For example, in astronomy an object can be observed with multiple filters, producing distinct time series of brightness that are not perfectly synchronized. To address this problem, a two-step procedure is proposed. The first step consists of pairing two time series by their observational times, in order to pair the measurements taken at the same time. In the case that for a given time there is a measurement in only one time series, a missing value is assigned in the remaining one. In the second step these missing values are imputed with an interpolation procedure.
To illustrate this procedure, we use two time series available in the iAR package. These time series, called cvnovag and cvnovar, correspond to the brightness magnitude of a cataclysmic variable object observed both in the g and r filters of the ZTF survey respectively and processed by the ALeRCE broker (Förster et al. 2021). Prior to pairing, a smooth trend is removed from each series using the loess function in order to obtain a stationary time series. This choice is not restrictive, and other smoothing or detrending methods available in R could be used instead (see, for instance, smooth.spline in base R, or functions available in packages such as mgcv (Wood (2001),Wood (2017)) and forecast (Hyndman and Khandakar 2008)). To pair these two time series, we use the method pairingits from the iAR package as follows:
This function causes that the datasets cvnovag and cvnovar of 67 and 65 observations respectively, to be merged into a single dataset with 97 observations which contains the measurements of both time series. Out of the total 97 observations in the merged dataset, 35 were matched based on their observation times, while in the remaining 62 observations there is a missing value in any of the two time series.
In order to fill these missing values we use the method interpolation. To implement this function we fit the BiAR model using the 35 observational times in which we have measures in the two time series as follows:
st=apply(pargr1[,c(1,4)],1,mean)
model_BiAR <- BiAR(times = st,series=pargr1[,c(2,5)],series_esd=pargr1[,c(3,6)])
model_BiAR@hessian=T
model_BiAR <- kalman(model_BiAR, niter = 20)
summary(model_BiAR,diagnostic=FALSE)BiAR model
Coefficients:
Estimate St. Error t value Pr(>|t|)
phiR 0.93 0.02 52.44 0.00
phiI 0.00 0.01 0.81 0.42
Information criteria:
AIC: -427.97
BIC: -424.86 [1] 0.9810805Note that the contemporary correlation between both time series estimated by the BiAR model is very high (rho=0.98). Using the estimated parameters by the BiAR model (phiR= 0.929 and phiI= 0.004), we estimate each missing value recursively, so that only one missing value is estimated per iteration. Note that by using the bivariate model to perform the smoothing, we achieve that the time series that has an observation at a given time helps the imputation of the time series that has the missing value at the same time. This is particularly useful for highly correlated series such as the one presented in this example.
Figures 11 a and 11 b show both time series after the pairing procedure, and Figures 11 c and 11 d show both time series after the smoothing procedure. Note that due to the high correlation between the raw time series, the imputed time series become very similar.
Figure 11: Figures a) and b) shows the time series cvnovag (green) and cvnovar (red) respectively. The vertical dotted lines show the observation times where a value must be imputed. The black dots are the estimated values at each of these observational times. Figures c) y d) shows both time series after the smoothing procedure.
7 Summary and discussion
A new R package called iAR is presented in this article. This package provides a flexible and consistent framework for modeling autoregressive processes with irregularly observed data. Its implementation using the S7 object-oriented system ensures a modular and extensible architecture, where each model class supports common methods for estimation, prediction, and interpolation. The use of S7 facilitates a coherent internal organization and simplifies future extensions to higher-dimensional models to be incorporated naturally.
This package implements novel models to fit irregularly observed time series. The methods implemented are the iAR, iAR-Gamma, iAR-T, CiAR and BiAR models. These models allow fitting a wide variety of autoregressive irregularly observed time series, some of which cannot be properly explained by models that assume regular or continuous time. For example, the CiAR model can fit irregularly observed time series with negative autocorrelation. Another feature is the ability to fit non-Gaussian data. The contexts in which each of these models is suitable are described throughout the article.
In the real-data applications presented in this paper, we showed that the models implemented in iAR are able to detect autoregressive dependence structures in irregularly observed time series and to generate interpolations in a model-based way. In this sense, the package provides an alternative to approaches that focus on interpolation of missing values, such as those implemented in imputeTS (Moritz and Bartz-Beielstein 2017), by explicitly modeling the underlying dependence structure rather than treating interpolation as a preprocessing step.
From a computational perspective, the current implementation is not parallelized, but the estimation procedures scale approximately linearly with the length of the time series. Empirical results show that the implemented functions allow estimating model parameters in less computational time than other R packages commonly used for this type of data.
The interface and functionality of iAR were developed with data analysts working with irregularly observed time series in mind, rather than exclusively with expert R programmers. Such data frequently arise in fields such as astronomy, remote sensing, climatology, finance, and medical studies, where observation times may be limited by external constraints. The proposed models are designed for stationary time series whose dependency structure is adequately described by the assumed autoregressive formulation, which is a natural limitation of the current framework.
Furthermore, an interactive application developed using the shiny package (Chang et al. 2020) allows users to fit the models implemented in iAR to their own data through a graphical interface. The application is publicly available at: https://iargem.shinyapps.io/iarFit/.
Overall, this package contributes novel methodological tools to the analysis of irregularly observed time series, a setting for which dedicated statistical software remains relatively scarce in the R ecosystem, despite its relevance across a wide range of applied disciplines (see, for instance, Ojeda et al. (2023b); Ojeda et al. (2023a); Godoy Pulecio and Ojeda Echeverri (2025)).
8 Computational details
The results in this paper were obtained using R~4.4.2 and the following packages: ggplot2~4.0.0, zoo~1.8-12, S7~0.2.0, Rcpp~1.1.0, Rdpack~2.6.2. R itself and all packages used are available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/.
Acknowledgments
The authors were supported by the ANID Millennium Science Initiative ICN12_009, awarded to the Millennium Institute of Astrophysics. F.E acknowledges the support of the ANID grant Fondecyt Iniciacion 11200590 and the CIRAS-AI Project code FIUF137139-USACH.
8.1 CRAN packages used
growth, yuima, carfima, sde, Rdrw, ctsem, lomb, RobPer, BINCOR, iAR, zoo, Rcpp, RcppArmadillo, mgcv, forecast, imputeTS, shiny, ggplot2, S7, Rdpack
8.2 CRAN Task Views implied by cited packages
Bayesian, ChemPhys, DifferentialEquations, Econometrics, Environmetrics, Finance, HighPerformanceComputing, MissingData, MixedModels, NetworkAnalysis, NumericalMathematics, Phylogenetics, Robust, Spatial, TeachingStatistics, TimeSeries, WebTechnologies