# Fitting univariate and multiviarate damped random walk processes 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. ```r 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. ## Returns 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

drw function