drw: Fitting Univariate and Multivariate Damped Random Walk...

In Rdrw: Univariate and Multivariate Damped Random Walk Processes

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)

Arguments

data1	An n_1 by 3 matrix or data frame for the first time series. The first column contains observation times, the second column contains measurements, and the third column contains measurement error standard deviations. If measurement error standard deviations are unavailable, use zeros in the third column.
data2	Optional second data set in the same format as data1.
data3	Optional third data set in the same format as data1.
data4	Optional fourth data set in the same format as data1.
data5	Optional fifth data set in the same format as data1.
data6	Optional sixth data set in the same format as data1.
data7	Optional seventh data set in the same format as data1.
data8	Optional eighth data set in the same format as data1.
data9	Optional ninth data set in the same format as data1.
data10	Optional tenth data set in the same format as data1. The current version allows up to ten time series data sets.
n.datasets	The number of time series data sets to be modeled jointly. This must be an integer from 1 to 10.
method	Either "mle" for maximum likelihood estimation or "bayes" for posterior sampling. Maximum likelihood estimation is computationally efficient and suitable for large-scale data analysis, while Bayesian posterior sampling is better suited to small-scale analyses requiring more careful uncertainty quantification.
poly.order	The order of the polynomial trend used for the mean function of each time series. The default is 1, corresponding to an intercept and a linear trend in normalized time. Use 0 for the constant-mean damped random walk model.
bayes.n.burn	The number of warm-up iterations when method = "bayes".
bayes.n.sample	The number of posterior samples saved after warm-up when method = "bayes".

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. The time variable u is the internally normalized observation time,

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.

Because the likelihood is evaluated on normalized time, the estimated tau values are also on the normalized-time scale. To convert them back to the original time scale, multiply by time.range = max(time) - min(time). For example, if the combined observation times span 2,000 days and tau = 0.05, then the corresponding timescale is 0.05 * 2000 = 100 days. The returned object includes this converted value as tau.days.

The polynomial coefficients in beta should be interpreted on the same normalized-time scale. 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.

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.

sigma2

Maximum likelihood estimates of the short-term variance parameters.

sigma

Maximum likelihood estimates of the short-term standard deviation parameters.

tau

Maximum likelihood estimates of the timescale parameters on the normalized-time scale.

tau.days

Maximum likelihood estimates of the timescale parameters on the original time scale. If the input observation times are measured in days, then tau.days is the estimated timescale in days. It is computed as tau * time.range.

rho

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

time.range

The range of the combined observation times, max(time) - min(time).

time.min

The minimum combined observation time.

data.comb

The merged data matrix used for fitting.

For method = "bayes", beta, sigma, tau, tau.days, 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.mle$tau.days


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)

Rdrw documentation built on Aug. 20, 2026, 5:09 p.m.