drw: Fitting univariate and multiviarate damped random walk processes

Description

The function drw fits univariate and multivariate damped random walk processes on multiple time series data sets possibly with known measurement error standard deviations via state-space representation. This function drw evaluates the resulting likelihood function of the model parameters via Kalman-filtering whose minimum complexity is linear in the number of unique observation times. The function returns the maximum likelihood estimates or posterior samples of the model parameters. For astronomical data analyses, users need to pay attention to loading the data because R's default is to load only seven effective digits; see details below.

Usage

drw(data1, data2, data3, data4, data5, 
    data6, data7, data8, data9, data10,
    n.datasets, method = "mle",
    bayes.n.burn, bayes.n.sample,
    mu.UNIFprior.range = c(-30, 30),
    tau.IGprior.shape = 1, tau.IGprior.scale = 1,
    sigma2.IGprior.shape = 1, sigma2.IGprior.scale = 1e-7)

Arguments

data1

An ($n_1$ by 3) matrix for time series data 1. The first column has $n_1$ observation times, the second column contains $n_1$ measurements (magnitudes), the third column includes $n_1$ measurement error standard deviations. If the measurement error standard deviations are unknown, the third column must be a vector of $n_1$ zeros.

data2

Optional if more than one time series data set are available. An ($n_2$ by 3) matrix for time series data 2. The first column has $n_2$ observation times, the second column contains $n_2$ measurements/observations/magnitudes, the third column includes $n_2$ measurement error standard deviations. If the measurement error standard deviations are unknown, the third column must be a vector of $n_2$ zeros.

data3

Optional if more than two time series data sets are available. An ($n_3$ by 3) matrix for time series data 3. The first column has $n_3$ observation times, the second column contains $n_3$ measurements/observations/magnitudes, the third column includes $n_3$ measurement error standard deviations. If the measurement error standard deviations are unknown, the third column must be a vector of $n_3$ zeros.

data4

Optional if more than three time series data sets are available. An ($n_4$ by 3) matrix for time series data 4. The first column has $n_4$ observation times, the second column contains $n_4$ measurements/observations/magnitudes, the third column includes $n_4$ measurement error standard deviations. If the measurement error standard deviations are unknown, the third column must be a vector of $n_4$ zeros.

data5

Optional if more than four time series data sets are available. An ($n_5$ by 3) matrix for time series data 5. The first column has $n_5$ observation times, the second column contains $n_5$ measurements/observations/magnitudes, the third column includes $n_5$ measurement error standard deviations. If the measurement error standard deviations are unknown, the third column must be a vector of $n_5$ zeros.

data6

Optional if more than five time series data sets are available. An ($n_6$ by 3) matrix for time series data 6. The first column has $n_6$ observation times, the second column contains $n_6$ measurements/observations/magnitudes, the third column includes $n_6$ measurement error standard deviations. If the measurement error standard deviations are unknown, the third column must be a vector of $n_6$ zeros.

data7

Optional if more than six time series data sets are available. An ($n_7$ by 3) matrix for time series data 7. The first column has $n_7$ observation times, the second column contains $n_7$ measurements/observations/magnitudes, the third column includes $n_7$ measurement error standard deviations. If the measurement error standard deviations are unknown, the third column must be a vector of $n_7$ zeros.

data8

Optional if more than seven time series data sets are available. An ($n_8$ by 3) matrix for time series data 8. The first column has $n_8$ observation times, the second column contains $n_8$ measurements/observations/magnitudes, the third column includes $n_8$ measurement error standard deviations. If the measurement error standard deviations are unknown, the third column must be a vector of $n_8$ zeros.

data9

Optional if more than eight time series data sets are available. An ($n_9$ by 3) matrix for time series data 9. The first column has $n_9$ observation times, the second column contains $n_9$ measurements/observations/magnitudes, the third column includes $n_9$ measurement error standard deviations. If the measurement error standard deviations are unknown, the third column must be a vector of $n_9$ zeros.

data10

Optional if more than nine time series data sets are available. An ($n_{10}$ by 3) matrix for time series data 10. The first column has $n_{10}$ observation times, the second column contains $n_{10}$ measurements/observations/magnitudes, the third column includes $n_{10}$ measurement error standard deviations. If the measurement error standard deviations are unknown, the third column must be a vector of $n_{10}$ zeros. The current version of package allows up to ten time series data sets. If users have more than ten time series data sets to be modeled simultaneously, please contact the authors.

n.datasets

The number of time series data sets to be modeled simultaneously. An integer value inclusively between 1 to 10. For example, if "n.datasets = 3", then users must enter data1, data2, and data3 as inputs of drw.

method

If method = "mle", it returns maximum likelihood estimates of the model parameters. If method = "bayes" it produces posterior samples of the model parameters.

bayes.n.burn

Required for method = "bayes". The number of warming-up iterations for a Markov chain Monte Carlo method.

bayes.n.sample

Required for method = "bayes". The size of a posterior sample for each parameter for a Markov chain Monte Carlo method.

mu.UNIFprior.range

Required for method = "bayes". The range of the Uniform prior on each long-term avaerage $\mu_j$of the process, where $j$ goes from 1 to the total number of time series data sets. The default range is ($-30, 30$) for astronomical applications.

tau.IGprior.shape

Required for method = "bayes". The shape parameter of the invserse-Gamma prior on each timescale $\tau_j$ of the process, where $j$ goes from 1 to the total number of time series data sets. The default shape parameter is one for astronomical applications.

tau.IGprior.scale

Required for method = "bayes". The scale parameter of the invserse-Gamma prior on each timescale $\tau_j$, where $j$ goes from 1 to the total number of time series data sets. The default scale parameter is one for astronomical applications.

sigma2.IGprior.shape

Required for method = "bayes". The shape parameter of the invserse-Gamma prior on each short-term variability (variance) $\sigma^2_j$, where $j$ goes from 1 to the total number of time series data sets. The default shape parameter is one for astronomical applications.

sigma2.IGprior.scale

Required for method = "bayes". The scale parameter of the invserse-Gamma prior on each short-term variability (variance) $\sigma^2_j$, where $j$ goes from 1 to the total number of time series data sets. The default shape parameter is 1e$-7$ for astronomical applications.

Value

The outcomes of drw are composed of:

mu: The maximum likelihood estimate(s) of the long-term average(s) if method is "mle", and the posterior sample(s) of the long-term average(s) if method is "bayes". In the former case (mle), it is a vector of length $k$, where $k$ is the number of time series data sets used. In the later case (bayes), it is an ($m$ by $k$) matrix where $m$ is the size of the posterior sample.
sigma: The maximum likelihood estimate(s) of the short-term variability (standard deviation) parameter(s) if method is "mle", and the posterior sample(s) of the short-term variability parameter(s) if method is "bayes". In the former case (mle), it is a vector of length $k$, where $k$ is the number of time series data sets used. In the later case (bayes), it is an ($m$ by $k$) matrix where $m$ is the size of the posterior sample.
tau: The maximum likelihood estimate(s) of the timescale(s) if method is "mle", and the posterior sample(s) of the timescale(s) if method is "bayes". In the former case (mle), it is a vector of length $k$, where $k$ is the number of time series data sets used. In the later case (bayes), it is an ($m$ by $k$) matrix where $m$ is the size of the posterior sample.
rho: Only when more than one time series data set are used. The maximum likelihood estimate(s) of the (cross-) correlation(s) if method is "mle", and the posterior sample(s) of the (cross-) correlation(s) if method is "bayes". In the former case (mle), it is a vector of length $k(k-1)/2$, where $k$ is the number of time series data sets used, i.e., $\rho_{12}, \ldots, \rho_{1k}, \rho_{23}, \ldots, \rho_{2k}, \ldots, \rho_{k-1, k}$. In the later case (bayes), it is an ($m$ by $k(k-1)/2$) matrix where $m$ is the size of the posterior sample.
mu.accept.rate: Only when method is "bayes". The MCMC acceptance rate(s) of the long-term average parameter(s).
sigma.accept.rate: Only when method is "bayes". The MCMC acceptance rate(s) of the short-term variability parameter(s).
tau.accept.rate: Only when method is "bayes". The MCMC acceptance rate(s) of the timescale(s).
rho.accept.rate: Only when more than one time series data set are used with method = "bayes". The MCMC acceptance rate(s) of the (cross-) correlation(s).
data.comb: The combined data set if more than one time series data set are used, and data1 if only one time series data set is used. This output is only available when method is set to "bayes".

Details

The multivariate damped random walk process ${\bf X}(t)$ is defined by the following stochastic differential equation: $$d {\bf X}(t) = -D_{{\bf\tau}}^{-1}({\bf X}(t) - {\bf\mu})dt + D_{{\bf\sigma}} d {\bf B}(t),$$ where ${\bf X}(t)=\{X_1(t), \ldots, X_k(t)\}$ is a vector of $k$ measurements/observations/magnitudes of the $k$ time series data sets in continuous time $t\in R$, $D_{\bf\tau}$ is a $k\times k$ diagonal matrix whose diagonal elements are $k$ timescales with each $\tau_j$ representing the timescale of the $j$-th time series data, ${\bf\mu}=\{\mu_1, \ldots, \mu_k\}$ is a vector for long-term averages of the $k$ time series data sets, $D_{{\bf\sigma}}$ is $k\times k$ diagonal matrix whose diagonal elements are short-term variabilities (standard deviation) of $k$ time series data sets, and finally ${\bf B}(t)=\{B_1(t), \ldots, B_k(t)\}$ is a vector for $k$ standard Brownian motions whose $k(k-1)/2$ pairwise correlations are modeled by correlation parameters $\rho_{jl}~(1\le j<l\le k)$ such that $dB_j(t)B_l(t) = \rho_{jl} dt$.

We evaluate this continuous-time process at $n$ discrete observation times ${\bf t} = \{t_1, \ldots, t_n\}$. The observed data ${\bf x} = \{x_1, \ldots, x_n\}$ are multiple time series data measured at irregularly spaced observation times ${\bf t}$ with possibly known measurement error standard deviations, $\delta=\{\delta_1, \ldots, \delta_n\}$. Since one or more time series observations can be measured at each observation time $t_i$, the length of a vector $x_i$ can be different, depending on how many time series observations are available at the $i$-th observation time. We assume that these observed data ${\bf x}$ are realizations of the latent time series data sets ${\bf X(t)} = \{{\bf X}(t_1), \ldots, {\bf X}(t_n)\}$ with Gaussian measurement errors whose variances are $\delta$. This is a typical setting of state-space modeling. We note that if the measurement error variances are unknown, $\delta$ must be set to zeros, which means that the observed data directly measure the latent values.

Please note that when astronomical time series data are loaded on R by read.table, read.csv, etc., some decimal places of the the observation times are automatically rounded because R's default is to load seven effective digits. For example, R will load the observation time 51075.412789 as 51075.41. This default will produce many ties in observation times even though there is actually no tie in observation times. To prevent this, please type "options(digits = 11)" before loading the data if the observation times are in seven effective digits.

References

Zhirui Hu and Hyungsuk Tak (2020+), "Modeling Stochastic Variability in Multi-Band Time Series Data," arXiv:2005.08049.

Examples

Run this code

# NOT RUN {
########## Fitting a univariate damped random walk process

##### Fitting a univariate damped random walk process based on a simulation

n1 <- 20 
# the number of observations in the data set

obs.time1 <- cumsum(rgamma(n1, shape = 3, rate = 1))
# the irregularly-spaced observation times 

y1 <- rnorm(n1)
# the measurements/observations/magnitudes

measure.error.SD1 <- rgamma(n1, shape = 0.01)
# optional measurement error standard deviations, 
# which is typically known in astronomical time series data
# if not known in other applications, set them to zeros, i.e.,
# measure.error.SD1 <- rep(0, n1)

data1 <- cbind(obs.time1, y1, measure.error.SD1)
# combine the single time series data set into an n by 3 matrix

# Note that when astronomical time series data are loaded on R (e.g., read.table, read.csv), 
# the digits of the observation times are typically rounded to seven effective digits.
# That means rounding may occur, which produces ties in observation times even though
# the original observation times are not the same.
# In this case, type the following code before loading the data.
# options(digits = 11)

res1.mle <- drw(data1 = data1, n.datasets = 1, method = "mle")
# obtain maximum likelihood estimates of the model parameters and 
# assign the result to object "res1.mle"

names(res1.mle)
# to see the maximum likelihood estimates,
# type "res1.mle$mu", "res1.mle$sigma", "res1.mle$tau"

# }
# NOT RUN {
res1.bayes <- drw(data1 = data1, n.datasets = 1, method = "bayes", 
                  bayes.n.burn = 10, bayes.n.sample = 10)
# }
# NOT RUN {
# obtain 10 posterior samples of each model parameter and
# save the result to object "res1.bayes"

# names(res1.bayes)
# to work on the posterior sample of each parameter, try 
# "res1.bayes$mu.accept.rate", "res1.bayes$sigma.accept.rate", "res1.bayes$tau.accept.rate"
# "hist(res1.bayes$mu)", "mean(res1.bayes$mu)", "sd(res1.bayes$mu)",
# "median(log(res1.bayes$sigma, base = 10))", 
# "quantile(log(res1.bayes$tau, base = 10), prob = c(0.025, 0.975))"


##### Fitting a multivariate damped random walk process based on simulations

n2 <- 10 
# the number of observations in the second data set

obs.time2 <- cumsum(rgamma(n2, shape = 3, rate = 1))
# the irregularly-spaced observation times of the second data set

y2 <- rnorm(n2)
# the measurements/observations/magnitudes of the second data set

measure.error.SD2 <- rgamma(n2, shape = 0.01)
# optional measurement error standard deviations of the second data set, 
# which is typically known in astronomical time series data
# if not known in other applications, set them to zeros, i.e.,
# measure.error.SD2 <- rep(0, n2)

data2 <- cbind(obs.time2, y2, measure.error.SD2)
# combine the single time series data set into an n by 3 matrix

# }
# NOT RUN {
res2.mle <- drw(data1 = data1, data2 = data2, n.datasets = 2, method = "mle")
# }
# NOT RUN {
# obtain maximum likelihood estimates of the model parameters and 
# assign the result to object "res2.mle"

# }
# NOT RUN {
res2.bayes <- drw(data1 = data1, data2 = data2, n.datasets = 2, method = "bayes", 
                  bayes.n.burn = 10, bayes.n.sample = 10)
# }
# NOT RUN {
# obtain 10 posterior samples of each model parameter and
# save the result to object "res2.bayes"

# names(res2.bayes)
# to work on the posterior sample of each parameter, try 
# "hist(res2.bayes$mu[, 1])", "colMeans(res2.bayes$mu)", "apply(res2.bayes$mu, 2, sd)",
# "hist(log(res2.bayes$sigma[, 2], base = 10))", 
# "apply(log(res2.bayes$sigma, base = 10), 2, median)", 
# "apply(log(res2.bayes$tau, base = 10), 2, quantile, prob = c(0.025, 0.975))"
# }

Run the code above in your browser using DataLab