Package {Rdrw}

Univariate and Multivariate Damped Random Walk Processes

Description

The Rdrw package provides tools for fitting and simulating univariate and multivariate damped random walk processes with optional measurement error standard deviations in a state-space framework. The damped random walk process is also known as an Ornstein–Uhlenbeck process or a first-order continuous-time autoregressive model, CAR(1) or CARMA(1, 0). The package uses Kalman filtering to evaluate the likelihood efficiently, with linear complexity in the number of unique observation times. Polynomial mean coefficients are handled conditionally within the Kalman recursion: they are profiled by generalized least squares for maximum likelihood estimation and sampled from their exact Gaussian conditional distribution for Bayesian inference. The main function drw fits the model and returns maximum likelihood estimates or posterior samples of the model parameters. The function drw.sim simulates time series from univariate or multivariate damped random walk processes.

Details

Author(s)

Zhirui Hu and Hyungsuk Tak

References

Zhirui Hu and Hyungsuk Tak (2020), "Modeling Stochastic Variability in Multi-Band Time Series Data," doi:10.3847/1538-3881/abc1e2.

Internal drw functions

Description

Internal drw functions

Details

These are not to be called by users.

Fitting Univariate and Multivariate Damped Random Walk Processes

Description

The function drw fits univariate and multivariate damped random walk processes to one or more irregularly spaced time series with possibly heteroscedastic measurement errors. The current implementation allows each time series to have a polynomial mean trend in normalized time. This feature is useful in applications such as stock price prediction, where a constant mean may be too restrictive. The likelihood is evaluated by Kalman filtering. Conditional on the damped random walk covariance parameters, the polynomial mean coefficients are linear Gaussian location parameters. For maximum likelihood estimation they are profiled analytically by generalized least squares within the Kalman recursion, and for Bayesian inference they are updated from their exact Gaussian conditional distribution. The function returns maximum likelihood estimates or posterior samples of the model parameters.

Usage

drw(data1, data2 = NULL, data3 = NULL, data4 = NULL, data5 = NULL,
    data6 = NULL, data7 = NULL, data8 = NULL, data9 = NULL,
    data10 = NULL, n.datasets, method = "mle", poly.order = 1,
    bayes.n.burn = 3000, bayes.n.sample = 1000,
    n.starts = 1, seed = NULL, optim.maxit = 3000,
    init.method = c("univariate", "legacy"))

Arguments

Details

The multivariate damped random walk process is also known as the Ornstein–Uhlenbeck process or a CARMA(1, 0) process. In the updated model, the long-term mean of the j-th time series is represented as

m_j(u) = beta_j0 + beta_j1 * u + ... + beta_jp * u^p,

where p is poly.order. For the polynomial mean only, the time coordinate u is defined as,

u = (t - t_min)/(t_max - t_min),

where t_min and t_max are the minimum and maximum combined observation times, respectively. Setting poly.order = 0 recovers the constant-mean specification in the original package.

The timescale tau is expressed in the same units as the supplied observation times. For example, with input times in days, tau = 100 means 100 days. The diffusion standard deviation sigma has units of the observations divided by the square root of the observation time unit.

The polynomial coefficients in beta are expressed in this dimensionless polynomial basis (unlike tau and sigma). For example, with poly.order = 1, the fitted mean is beta_j0 + beta_j1 * u, so beta_j1 is the total linear change over the normalized interval [0, 1], not the daily slope.

For computational efficiency, the polynomial coefficients are not included in the nonlinear optimizer used for maximum likelihood estimation. At each proposed set of covariance parameters, the Kalman recursion accumulates the Gaussian sufficient statistics for the polynomial coefficients and solves for their conditional generalized least squares estimates exactly. The coefficients remain estimated model parameters and are returned in the same order as in previous versions.

For Bayesian inference, the same Kalman sufficient statistics define the exact Gaussian conditional distribution of the polynomial coefficients under the package's diffuse Gaussian prior. The sampler draws these coefficients directly and caches the current Kalman sufficient statistics and log likelihood during the Metropolis–Hastings updates of the covariance parameters, avoiding redundant likelihood evaluations.

Initialization and optimization. By default, each marginal series is fitted as a univariate OU/DRW model before joint inference. Each preliminary fit uses three observation-adaptive timescale initializations spanning representative within-band cadence and observational baseline; initial diffusion variance is matched to the measurement-error-adjusted empirical stationary variance. The preliminary fits provide only starting values, not fixed parameters, likelihood penalties, or Bayesian priors. If a univariate fit is unavailable, a cadence/baseline-based fallback is used.

For joint MLE, n.starts = 1 uses the univariate warm starts and an equicorrelation matrix with all pairwise correlations equal to 0.5. With two or more starts, the second uses independence; the third, when requested, uses a shrunken correlation estimate of overlapping standardized marginal Kalman innovations (or a random fallback). Subsequent starts perturb the marginal fits and explore observation-adaptive random timescales and valid random correlation matrices. The innovation-based correlations are numerical initializations, not estimates of the driving-noise correlations. Correlations are parameterized using Cholesky partial correlations, guaranteeing positive definiteness. All starts optimize the original joint likelihood and the best solution is refined. init.method = "legacy" skips marginal fits and uses a fixed normalized-time initialization. n.starts controls joint MLE only.

For method = "bayes", the same marginal warm starts and valid equicorrelation initialization seed the MCMC directly, without requiring a preliminary joint MLE. The correlation sampler uses bounded unconstrained Cholesky partial correlations, with a uniform prior on these transformed coordinates (not a uniform prior on the pairwise correlations). Symmetric random-walk proposals preserve the target posterior on this parameterization. Other covariance parameters retain bounded uniform priors on their logarithmic coordinates, and the polynomial mean coefficients retain a diffuse Gaussian prior. Data-informed initial states do not enter the posterior prior or likelihood. Multiple chains and convergence checks are recommended.

The result omits the merged data matrix; marginal warm-start diagnostics are in initialization. Estimated timescales are reported in the observation time units. The diffusion standard deviation sigma is distinct from the stationary standard deviation: the latter equals sigma * sqrt(tau / 2).

For multivariate fitting, the individual data sets are merged by observation time. Missing measurements at a given time are allowed and are handled by the Kalman filter. The input measurement-error column should contain standard deviations, not variances.

When astronomical time series data are loaded into R by read.table, read.csv, or related functions, decimal places of observation times may be rounded because of R's display and input defaults. If high time precision is needed, consider setting options(digits = 11) before loading the data.

Value

The returned object is a list. For method = "mle", it contains:

beta

Maximum likelihood estimates of the polynomial trend coefficients. The coefficients are returned as a vector ordered by time series. For each time series, the coefficients correspond to powers 0 through poly.order of normalized time, u = (t - t_min)/(t_max - t_min). Thus, beta_j0 is the fitted mean at the first normalized time point, and higher-order coefficients describe changes over the unit interval of normalized time, not per day.

sigma

Estimated diffusion standard deviations, in units of the observations divided by the square root of the observation time unit.

tau

Estimated characteristic timescales, in the same units as the observation times.

rho

Maximum likelihood estimates of the cross-correlation parameters when n.datasets > 1; otherwise numeric(0).

For method = "bayes", beta, sigma, tau, and rho contain posterior samples, with one row per saved MCMC iteration. The output also includes sigma.accept.rate, tau.accept.rate, and rho.accept.rate.

Author(s)

Zhirui Hu and Hyungsuk Tak

References

Zhirui Hu and Hyungsuk Tak (2020), "Modeling Stochastic Variability in Multi-Band Time Series Data," doi:10.3847/1538-3881/abc1e2.

Examples

set.seed(1)
n1 <- 20
obs.time1 <- cumsum(rgamma(n1, shape = 3, rate = 1))
y1 <- rnorm(n1)
measure.error.SD1 <- rgamma(n1, shape = 0.01)
data1 <- cbind(obs.time1, y1, measure.error.SD1)

res1.mle <- drw(data1 = data1, n.datasets = 1, method = "mle", poly.order = 1)
names(res1.mle)
res1.mle$beta
res1.mle$sigma
res1.mle$tau


res1.bayes <- drw(data1 = data1, n.datasets = 1, method = "bayes",
                  poly.order = 1, bayes.n.burn = 10, bayes.n.sample = 10)


n2 <- 15
obs.time2 <- cumsum(rgamma(n2, shape = 3, rate = 1))
y2 <- rnorm(n2)
measure.error.SD2 <- rgamma(n2, shape = 0.01)
data2 <- cbind(obs.time2, y2, measure.error.SD2)


res2.mle <- drw(data1 = data1, data2 = data2, n.datasets = 2,
                method = "mle", poly.order = 1, n.starts = 5)

Simulating Univariate and Multivariate Damped Random Walk Processes

Description

The function drw.sim simulates one or more time series from univariate or multivariate damped random walk processes.

Usage

drw.sim(time, n.datasets, measure.error.SD, mu, sigma, tau, rho)

Arguments

Details

Given observation times and model parameter values mu, sigma, tau, and optionally rho, this function simulates k time series, optionally including known measurement error standard deviations.

Value

The output of drw.sim is a list containing:

x

An n by k matrix of simulated time series. Each column corresponds to one simulated time series.

Author(s)

Zhirui Hu and Hyungsuk Tak

References

Zhirui Hu and Hyungsuk Tak (2020), "Modeling Stochastic Variability in Multi-Band Time Series Data," doi:10.3847/1538-3881/abc1e2.

Examples

########## Simulating a multivariate damped random walk process

n <- 100
k <- 5
obs.time <- cumsum(rgamma(n, shape = 3, rate = 1))

tau <- 100 + 20 * (1:5)
sigma <- 0.01 * (1:5)
mu <- 17 + 0.5 * (1:5)

rho.m <- matrix(0, k, k)
for (i in 1:k) {
  for (j in 1:k) {
    rho.m[i, j] <- 1.1^(-abs(i - j))
  }
}

rho <- rho.m[upper.tri(rho.m)]

measure.error.band <- c(0.010, 0.014, 0.018, 0.022, 0.026)
measure.error <- NULL
for (i in 1:k) {
  measure.error <- cbind(measure.error, rnorm(n, measure.error.band[i], 0.002))
}

x <- drw.sim(time = obs.time, n.datasets = 5, measure.error.SD = measure.error,
             mu = mu, sigma = sigma, tau = tau, rho = rho)

plot(obs.time, x[, 1], xlim = c(min(obs.time), max(obs.time)), ylim = c(17, 20),
     xlab = "time", ylab = "observation")
points(obs.time, x[, 2], col = 2, pch = 2)
points(obs.time, x[, 3], col = 3, pch = 3)
points(obs.time, x[, 4], col = 4, pch = 4)
points(obs.time, x[, 5], col = 5, pch = 5)

########## Simulating a univariate damped random walk process

x <- drw.sim(time = obs.time, n.datasets = 1, measure.error.SD = measure.error[, 1],
             mu = mu[1], sigma = sigma[1], tau = tau[1])
plot(obs.time, x)