# Stochastic Volatility Regression for
Functional Data Dynamics

Bin Zhu and David B. Dunson∗

††footnotetext: ∗Bin Zhu is Tenure-Track Principal Investigator, Biostatistics
Branch, Division of Cancer Epidemiology and Genetics, National Cancer
Institute, National Institutes of Health, Rockville, MD 20852 (Email:
_bin.zhu@nih.gov_). David B. Dunson is Professor, Department of Statistical
Science, Duke University, Durham, NC 27708, (Email: _dunson@stat.duke.edu_).

Abstract

Although there are many methods for functional data analysis (FDA), little
emphasis is put on characterizing variability among volatilities of individual
functions. In particular, certain individuals exhibit erratic swings in their
trajectory while other individuals have more stable trajectories. There is
evidence of such volatility heterogeneity in blood pressure trajectories
during pregnancy, for example, and reason to suspect that volatility is a
biologically important feature. Most FDA models implicitly assume similar or
identical smoothness of the individual functions, and hence can lead to
misleading inferences on volatility and an inadequate representation of the
functions. We propose a novel class of FDA models characterized using
hierarchical stochastic differential equations. We model the derivatives of a
mean function and deviation functions using Gaussian processes, while also
allowing covariate dependence including on the volatilities of the deviation
functions. Following a Bayesian approach to inference, a Markov chain Monte
Carlo algorithm is used for posterior computation. The methods are tested on
simulated data and applied to blood pressure trajectories during pregnancy.
Key words: Diffusion process; Gaussian process; State space model; Stochastic
differential equation; Stochastic dynamic model; Stochastic functional data
analysis.

##  1 Introduction

Multi-subject functional data arise frequently in many fields of research,
including epidemiology, clinical trials and environmental health. Sequential
observations are collected over time for multiple subjects, and can be treated
as being generated from a smooth trajectory contaminated with noise. There are
a rich variety of methods available for characterizing variability and
covariate dependence in functional data ranging from hierarchical basis
expansions to functional principal components analysis (FPCA). In defining
models for functional data and in interpreting variability in trajectories,
either unexplained or due to covariates, the emphasis has been on differences
in the level and local trends. Dynamic features, such as velocity,
acceleration and especially volatility, are also important but receive much
less attention.

Analysis of functional data dynamics studies temporal changes in trajectories
and effects of covariates on these changes. For example, Wang et al. (2008)
used linear differential equations to model price velocity and acceleration in
eBay auctions and explored the auction subpopulation effect. Müller and Yao
(2010) modeled the velocity of online auction bids using empirical stochastic
differential equations with time-varying coefficients and a smooth drift
process. Zhu et al. (2011) inferred the rate functions of prostate-specific
antigen profiles using the proposed semiparametric stochastic velocity model,
in which rate functions are regarded as realizations of Ornstein-Uhlenbeck
processes dependent on covariates of interest.

This article investigates a different dynamic feature, the volatility, which
measures the conditional variance of trajectory changes over an infinitesimal
time interval. We propose a stochastic volatility regression (SVR) model with
Gaussian process (GP) priors used for the group mean and subject specific
deviation functions through stochastic differential equations (SDEs). We
further accommodate inference on covariate effects on volatility through
allowing the diffusion term of SDEs for deviation functions to depend on
covariates. Although volatility has been extensively studied through
stochastic volatility (SV) models in finance (Heston, 1993; Jacquier et al.,
2002; Shephard, 2005; Barndorff-Nielsen and Shephard, 2012), the setting,
model specifications and data features are distinct from ours. SV models
typically focus on a single volatility process which is time-varying and
associated with a price process for high-frequency finance data. More relevant
is the literature on multivariate SV models; for recent references, refer to
Loddo et al. (2011), Van Es and Spreij (2011), Ishihara and Omori (2012) and
Durante et al. (2012).

This setting differs from ours in that the focus is on multivariate time
series modeling instead of functional data analysis, with interest in the
joint volatility dynamics over time for the different assets. In contrast, we
are interested in studying variation across individuals in a time-constant
subject-specific volatility; that is, certain subjects may have very smooth
trajectories while other subjects have erratic trajectories. It is our
conjecture that such volatility heterogeneity is common in biomedical
settings, but is overlooked in analyzing data with models that implicitly
prescribe a single level of smoothness for all subjects. As data are sparse
and irregularly spaced in most studies, it is not surprising such behavior is
overlooked. However, the volatility in a biomarker may be as important or more
important than the overall level and trend in the biomarker. We provide
motivation through the following longitudinal blood pressure data set.

The Healthy Pregnancy, Healthy Baby Study (HPHB, Miranda et al., 2009)
collected longitudinal blood pressure (BP) measurements for pregnant women.
Blood pressures are measured at irregularly spaced times during the second and
third trimesters with the number of measurements per subject varying from 9 to
19. We are interested in estimating subject-specific volatilities of BP
trajectories and in identifying covariates associated with the volatility.
Figure 1LABEL:sub@fig:1a plots mean arterial pressure (MAP) trajectories for
twenty randomly selected normal women and women with preeclampsia,
respectively. Clearly the MAP trajectories among the preeclampsia group are
more wiggly than the ones in the normal group, which is also implied by
boxplots of log-transformed empirical volatilities in Figure
1LABEL:sub@fig:1b. To explore volatility differences among various groups in
addition to preeclampsia, we apply normal linear regression for log-
transformed empirical volatilities with the covariates race, mother’s age,
obesity, preeclampsia, parity and smoking. The results suggest that
preeclampsia and smoking (p-values 0.0005 and 0.002) are associated with
empirical volatility. This is a two-stage approach, which is useful as an
exploratory tool but ignores measurement errors and uncertainty in volatility
estimation.

(a)

(b)

Figure 1: (a) Mean arterial pressure (MAP) trajectories for twenty randomly
selected normal women and women with preeclampsia; (b) Log-transformed
empirical volatilities for women in the normal group and preeclampsia group.
$Y_{ij}$ denotes blood pressure for the $i$th woman at time $t_{ij}$, and
$U_{ij}=Y_{ij}-\bar{Y}_{j}$ is the deviation from the group mean blood
pressure $\bar{Y}_{j}$. The empirical volatility measures the fluctuation of
trajectories empirically and is defined as
$\frac{1}{n_{i}}\sum_{j=1}^{n_{i}-1}\frac{(U_{i,j+1}-U_{i,j})^{2}}{t_{i,j+1}-t_{i,j}}$
with $n_{i}$ the number of observations for the $i$th woman.

Additionally, empirical volatilities in Figure 1LABEL:sub@fig:1b are
heterogeneous even within the normal or preeclampsia group. This heterogeneity
will be largely omitted when we apply FDA methods with identical or similar
smoothness for individual functions within a group. Consequently, the wiggly
trajectories will be over-smoothed while the smooth trajectories will be
under-smoothed. We can potentially estimate the individual trajectories
separately but it is well known that borrowing of information will
dramatically improve performance for sparse functional data. In addition,
separate estimation does not allow for inferences on covariate effects and
unexplained variability in volatility.

As for the clinical question addressed, the previous FDA methods mainly focus
on the shift of blood pressure level and ignore examining the volatility of
blood pressure, which measures the haemodynamic stability and is crucial for
cardiovascular health. For example, a recent study shows that blood pressure
stability rather than blood pressure level is associated with increased
survival among patients on hemodialysis (Raimann et al., 2012). For the HPHB
study, we observe that preeclampsia is commonly accompanied by blood pressure
over-swinging. The joint effect of high blood pressure level and large
volatility may lead to the adverse birth outcomes, such as low birth weight
and preterm birth.

The remainder of the article is organized as follows. Section 2 specifies the
SVR model and discusses its properties. Section 3 develops an efficient Markov
chain Monte Carlo algorithm for posterior inference. Section 4 presents
simulation studies and the proposed method is applied to a real dataset in
Section 5. Finally, section 6 contains concluding remarks and future possible
extensions.

##  2 Stochastic Volatility Regression Model

###  2.1 The Model Specification

Suppose that $Y_{i}(t)$, $i=1,2,\dots,m$, is the observation of the $i$th
subject at time
$t\in\mathcal{T}_{i}=\\{t_{i,1},t_{i,2},\cdots,t_{i,n_{i}}<t_{U}\\}$ with
$\mathcal{T}_{i}$ the set of observation times before time $t_{U}$ for the
$i$th subject. We specify an observation equation for $Y_{i}(t)$ as

$Y_{i}(t)=M_{k_{i}}(t)+U_{i}(t)+\varepsilon_{i}(t),$ (1)

where $Y_{i}(t)$ is contaminated by the measurement error $\varepsilon_{i}(t)$
following a one-dimensional normal distribution with mean $0$ and variance
$\sigma^{2}_{\varepsilon}$. Assuming the $i$th subject belongs to the
${k_{i}}$th group (e.g. by race or treatment) with
${k_{i}}\in\\{1,2,\cdots,g\\}$, we include a ${k_{i}}$th group mean function
$M_{k_{i}}(t)=\textsf{E}\\{Y_{i}(t)\mid M_{k_{i}}(t)\\}$ in the observation
equation. In addition, the trajectory of the $i$th subject will unlikely
coincide with $M_{k_{i}}(t)$ and therefore the departure from $M_{k_{i}}(t)$
is addressed and represented by the subject-specific deviation function
$U_{i}(t)$ with $\textsf{E}\\{U_{i}(t)\\}=0$.

The volatility of the $i$th subject is defined as the conditional variance of
the $(q-1)$th order derivative of $U_{i}(t)$ over an infinitesimal time
interval. Namely, we denote the volatility
$\sigma^{2}_{U_{i}}=\mathop{\lim}\limits_{h\to
0}h^{-1}\textsf{E}\left[\left\\{D^{q-1}U_{i}(t+h)-D^{q-1}U_{i}(t)\right\\}^{2}\mid
D^{q-1}U_{i}(t)\right]$ with differential operator
$D^{q}=\frac{d^{q}}{dt^{q}}$. As volatility approaches zero, $U_{i}(t)$ would
be a roughly flat line. In contrast, increasing the value of volatility would
lead to a more wiggly $U_{i}(t)$ with a larger magnitude of fluctuation around
$M_{k_{i}}(t)$.

We specify Gaussian process priors for $M_{k_{i}}(t)$ and $U_{i}(t)$ using
SDEs which incorporate the group and individual volatilities
$\sigma^{2}_{M_{k_{i}}}$ and $\sigma^{2}_{U_{i}}$ :

$\displaystyle D^{p}M_{k_{i}}(t)$
$\displaystyle=\sigma_{M_{k_{i}}}\dot{W}_{k_{i}}(t),$ (2) $\displaystyle
D^{q}U_{i}(t)$ $\displaystyle=\sigma_{U_{i}}\dot{W}^{\prime}_{i}(t),$ (3)

where $p,q\in\mathbb{N}\geq 1$ and
$\sigma_{M_{k_{i}}},\sigma_{U_{i}}\in\mathbb{R}^{+}$; $\dot{W}_{k_{i}}(t)$ and
$\dot{W}^{\prime}_{i}(t)$ are independent Gaussian white noise processes with
$\textsf{E}\\{\dot{W}_{k_{i}}(t)\\}=\textsf{E}\\{\dot{W}^{\prime}_{i}(t)\\}=0$
and covariance function
$\textsf{E}\\{\dot{W}_{k_{i}}(t)\dot{W}_{k_{i}}(t^{\prime})\\}=\textsf{E}\\{\dot{W}^{\prime}_{i}(t)\dot{W}^{\prime}_{i}(t^{\prime})\\}=\delta(t-t^{\prime})$,
a delta function. We denote
$\boldsymbol{M}_{{k_{i}}0}=\\{M_{k_{i}}(0),D^{1}M_{k_{i}}(0),\dots,D^{p-1}M_{k_{i}}(0)\\}$
and $\boldsymbol{U}_{i0}=\\{U_{i}(0),D^{1}U_{i}(0),\dots,D^{q-1}U_{i}(0)\\}$
as the initial values of $M_{k_{i}}(t)$ and $U_{i}(t)$ as well as their
derivatives till orders $q-1$ and $p-1$ respectively. The volatility
$\sigma^{2}_{U_{i}}$ in SDE (3) is allowed to vary between subjects and across
covariates. In this article, we focus on a simple transformed mean
relationship, namely
$\log(\sigma^{2}_{U_{i}})\sim\mathsf{N}_{1}(\boldsymbol{x}^{\prime}_{i}\boldsymbol{\beta},\sigma^{2})$,
which can be extended to the more complex specifications with less restrictive
assumptions and to high-dimensional covariates.

The mean and covariance functions of Gaussian process priors for
$M_{k_{i}}(t)$ and $U_{i}(t)$ can be obtained by applying stochastic
integration to SDEs (2) and (3), resulting in the following lemma.

###### Lemma 1

$M_{k}(t)$, $k=1,2,\dots,g$, and $U_{i}(t)$, $i=1,2,\dots,n$, are the
summations of mutually independent Gaussian processes written as
$M_{k}(t)=M_{k0}(t)+M_{k1}(t)$ and $U_{i}(t)=U_{i0}(t)+U_{i1}(t)$ with
corresponding mean functions
$\textsf{E}\left\\{M_{k0}(t)\right\\}=\textsf{E}\left\\{M_{k1}(t)\right\\}=\textsf{E}\left\\{U_{i0}(t)\right\\}=\textsf{E}\left\\{U_{i1}(t)\right\\}=0$
and covariance functions

$\displaystyle\mathcal{K}_{M_{k0}}(s,t)$
$\displaystyle=\sigma^{2}_{M_{0}}\mathcal{R}_{M_{0}}(s,t)=\sigma^{2}_{M_{0}}\sum_{l=0}^{q-1}\phi_{l}(s)\phi_{l}(t),$
$\displaystyle\mathcal{K}_{M_{k1}}(s,t)$
$\displaystyle=\sigma^{2}_{M_{k}}\mathcal{R}_{M_{1}}(s,t)=\sigma^{2}_{M_{k}}\int_{\mathcal{T}}G_{q}(s,u)G_{q}(t,u)du,$
$\displaystyle\mathcal{K}_{U_{i0}}(s,t)$
$\displaystyle=\sigma^{2}_{U_{0}}\mathcal{R}_{U_{0}}(s,t)=\sigma^{2}_{U_{0}}\sum_{l=0}^{p-1}\phi_{l}(s)\phi_{l}(t),$
$\displaystyle\mathcal{K}_{U_{i1}}(s,t)$
$\displaystyle=\sigma^{2}_{U_{i}}\mathcal{R}_{U_{1}}(s,t)=\sigma^{2}_{U_{i}}\int_{\mathcal{T}}G_{p}(s,u)G_{p}(t,u)du,$

respectively, where $\phi_{l}(t)=\frac{t^{l}}{l!}$,
$G_{q}(s,u)=\frac{(s-u)_{+}^{q-1}}{(q-1)!}$ and
$s,t,u\in\mathcal{T}=[0,t_{U}]$.

Hence, we can represent the prior of $M_{k_{i}}(t)+U_{i}(t)$ as a hierarchical
Gaussian process,

$\displaystyle M_{k_{i}}(t)+U_{i}(t)\mid M_{k_{i}}(t)$
$\displaystyle\sim\mathcal{GP}(M_{k_{i}}(t),\mathcal{K}_{U_{i0}}(s,t)+\mathcal{K}_{U_{i1}}(s,t)),$
$\displaystyle M_{k_{i}}(t)$
$\displaystyle\sim\mathcal{GP}(0,\mathcal{K}_{M_{k0}}(s,t)+\mathcal{K}_{M_{k1}}(s,t)),$

where $\mathcal{GP}(M(t),\mathcal{K}(s,t))$ denotes a Gaussian process with
mean function $M(t)$ and covariance function $\mathcal{K}(s,t)$. Different
from the previous hierarchical Gaussian process prior (Park and Choi, 2010),
in which the covariance function is modeled as a squared exponential kernel
and is identical across the subjects within a group, here
$\mathcal{K}_{U_{i0}}(s,t)+\mathcal{K}_{U_{i1}}(s,t)$ is subject specific and
depends on covariates through $\sigma^{2}_{U_{i}}$.

To carry out Bayesian inference, we further specify the following prior
distributions for the hyperparameters. In particular,
$\boldsymbol{M}_{{k_{i}}0}\sim\mathsf{N}_{p}(\boldsymbol{0},\sigma^{2}_{M_{0}}\boldsymbol{I})$
with $\sigma^{2}_{M_{0}}=10^{4}$,
$\boldsymbol{U}_{i0}\sim\mathsf{N}_{q}(\boldsymbol{0},\sigma^{2}_{U_{0}}\boldsymbol{I})$,
$\sigma^{2}_{\varepsilon}\sim\textsf{invGamma}(a,b)$,
$\sigma^{2}_{M_{k}}\sim\textsf{invGamma}(a,b)$ and
$\sigma^{2}_{U_{0}}\sim\textsf{invGamma}(a,b)$, where $\textsf{invGamma}(a,b)$
denotes the inverse gamma distribution with shape parameter $a$ and scale
parameter $b$. The $\boldsymbol{\beta}$ and $\sigma^{2}$ follow the
independent Jeffreys’ prior, $f(\boldsymbol{\beta},\sigma^{2})\propto
1/\sigma^{2}$.

###  2.2 Double-Penalized Smoothing Spline

It is well known that the Bayes estimate with the integrated Wiener process
prior is identical to the smoothing spline estimate (Wahba, 1990). By similar
arguments, we can show that when the volatilities are given and
$\sigma^{2}_{M_{0}}$ and $\sigma^{2}_{U_{0}}$ go to infinity, the posterior
means of $M_{k}(t)$ and $U_{i}(t)$ are equivalent to the double-penalized
smoothing spline $\hat{M}_{k}(t)+\hat{U}_{i}(t)$, which is the minimizer of
the double-penalized sum-of-squares,

$\displaystyle\textsf{DPSS}=$
$\displaystyle\sum_{i=1}^{m}\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\left\\{Y(t_{ij})-M_{k_{i}}(t_{ij})-U_{i}(t_{ij})\right\\}^{2}+$
(4)
$\displaystyle\sum_{k=1}^{g}\lambda_{M_{k}}\int_{\mathcal{T}}\left\\{D^{p}M_{k}(t)\right\\}^{2}dt+\sum_{i=1}^{m}\lambda_{U_{i}}\int_{\mathcal{T}}\left\\{D^{q}U_{i}(t)\right\\}^{2}dt,$

where penalty terms $\int_{\mathcal{T}}\left\\{D^{p}M_{k}(t)\right\\}^{2}dt$
and $\int_{\mathcal{T}}\left\\{D^{q}U_{i}(t)\right\\}^{2}dt$ penalize the
roughness of $M_{k}(t)$ and $U_{i}(t)$ respectively, where the smoothness and
the fidelity to data are balanced by the smoothing parameters
$\lambda_{M_{k}}=\displaystyle\sum_{i:k_{i}=k}\frac{\sigma^{2}_{\varepsilon}}{n_{i}\sigma^{2}_{M_{k}}}$
and
$\lambda_{U_{i}}=\frac{\sigma^{2}_{\varepsilon}}{n_{i}\sigma^{2}_{U_{i}}}$.
Expression for $\hat{M}_{k}(t)$ and $\hat{U}_{i}(t)$ can be obtained
explicitly, as shown in the following Theorem.

###### Theorem 1

The smoothing splines $\hat{M}_{k}(t)$ and $\hat{U}_{i}(t)$ with
$t\in\mathcal{T}$ minimize the double-penalized sum-of-squares (4) and have
the forms

$\displaystyle\hat{M}_{k}(t)$
$\displaystyle=\sum_{l=0}^{p-1}\mu_{kl}\phi_{l}(t)+\sum_{j=1}^{n}\nu_{kj}\mathcal{R}_{M_{1}}(t_{j},t)=\boldsymbol{\mu}^{\prime}_{k}\boldsymbol{\phi}_{\mu}(t)+\boldsymbol{\nu}^{\prime}_{k}\boldsymbol{R}_{M_{1}}(t)$
$\displaystyle\hat{U}_{i}(t)$
$\displaystyle=\sum_{l=0}^{q-1}\alpha_{il}\phi_{l}(t)+\sum_{j=1}^{n_{i}}\gamma_{ij}\mathcal{R}_{U_{1}}(t_{ij},t)=\boldsymbol{\alpha}^{\prime}_{i}\boldsymbol{\phi}_{\alpha}(t)+\boldsymbol{\gamma}^{\prime}_{i}\boldsymbol{R}_{U_{i1}}(t)$

where
$\boldsymbol{\mu}_{k}=\\{\mu_{k0},\mu_{k1},\cdots,\mu_{k(p-1)}\\}^{\prime}$,
$\boldsymbol{\nu}_{k}=(\nu_{k1},\nu_{k2},\cdots,\nu_{kn})^{\prime}$,
$\boldsymbol{\alpha}_{i}=\left\\{\alpha_{i0},\alpha_{i1},\cdots,\alpha_{i(q-1)}\right\\}^{\prime}$
and
$\boldsymbol{\gamma}_{i}=(\gamma_{i1},\gamma_{i2},\cdots,\gamma_{in_{i}})^{\prime}$
are the coefficients for the bases

$\displaystyle\boldsymbol{\phi}_{\mu}(t)$
$\displaystyle=\\{\phi_{0}(t),\phi_{1}(t),\cdots,\phi_{p-1}(t)\\}^{\prime},\quad\boldsymbol{R}_{M_{1}}(t)=\\{\mathcal{R}_{M_{1}}(t_{1},t),\mathcal{R}_{M_{1}}(t_{2},t),\cdots,\mathcal{R}_{M_{1}}(t_{n},t)\\}^{\prime},$
$\displaystyle\boldsymbol{\phi}_{\alpha}(t)$
$\displaystyle=\\{\phi_{0}(t),\phi_{1}(t),\cdots,\phi_{q-1}(t)\\}^{\prime},\quad\boldsymbol{R}_{U_{i1}}(t)=\\{\mathcal{R}_{U_{1}}(t_{i1},t),\mathcal{R}_{U_{1}}(t_{i2},t),\cdots,\mathcal{R}_{U_{1}}(t_{in_{i}},t)\\}^{\prime},$

with
$t_{j}\in\mathcal{T}_{m}=\cup_{i=1}^{m}\mathcal{T}_{i}=\\{t_{j},\;j=1,2,\dots,n\\}$,
the index set of unique observation times among all $m$ subjects.

Given $\hat{M}_{k}(t)$ and $\hat{U}_{i}(t)$, the double-penalized sum-of-
squares (4) can be written as

$\displaystyle\textsf{DPSS}=$
$\displaystyle\sum_{i=1}^{m}\frac{1}{n_{i}}(\boldsymbol{Y}_{i}-\boldsymbol{\Delta}_{i}\boldsymbol{\phi}_{\mu}\boldsymbol{\mu}_{k_{i}}-\boldsymbol{\Delta}_{i}\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k_{i}}-\boldsymbol{\phi}_{\alpha_{i}}\boldsymbol{\alpha}_{i}-\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i})^{\prime}\times$
(5)
$\displaystyle(\boldsymbol{Y}_{i}-\boldsymbol{\Delta}_{i}\boldsymbol{\phi}_{\mu}\boldsymbol{\mu}_{k_{i}}-\boldsymbol{\Delta}_{i}\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k_{i}}-\boldsymbol{\phi}_{\alpha_{i}}\boldsymbol{\alpha}_{i}-\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i})+$
$\displaystyle\sum_{k=1}^{g}\lambda_{M_{k}}\boldsymbol{\nu}^{\prime}_{k}\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k}+\sum_{i=1}^{m}\lambda_{U_{i}}\boldsymbol{\gamma}^{\prime}_{i}\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i},$

where

$\displaystyle\boldsymbol{Y}_{i}$
$\displaystyle=\\{Y(t_{i1}),Y(t_{i2}),\cdots,Y(t_{in_{i}})\\}^{\prime},\quad\boldsymbol{\Delta}_{i}=(\delta_{jj^{\prime}})_{n_{i}\times
n},$ $\displaystyle\boldsymbol{\phi}_{\mu}$
$\displaystyle=\\{\boldsymbol{\phi}_{\mu}(t_{1}),\boldsymbol{\phi}_{\mu}(t_{2}),\cdots,\boldsymbol{\phi}_{\mu}(t_{n})\\}^{\prime},\quad\boldsymbol{R}_{M_{1}}=\\{\boldsymbol{R}_{M_{1}}(t_{1}),\boldsymbol{R}_{M_{1}}(t_{2}),\cdots,\boldsymbol{R}_{M_{1}}(t_{n})\\},$
$\displaystyle\boldsymbol{\phi}_{\alpha_{i}}$
$\displaystyle=\\{\boldsymbol{\phi}_{\alpha}(t_{i1}),\boldsymbol{\phi}_{\alpha}(t_{i2}),\cdots,\boldsymbol{\phi}_{\alpha}(t_{in_{i}})\\}^{\prime},\quad\boldsymbol{R}_{U_{i1}}=\\{\boldsymbol{R}_{U_{i1}}(t_{i1}),\boldsymbol{R}_{U_{i1}}(t_{i2}),\cdots,\boldsymbol{R}_{U_{i1}}(t_{in_{i}})\\}$

with $\delta_{jj^{\prime}}=1$ if $i$th subject has an observation at time
$t_{ij}=t_{j^{\prime}}$, $t_{ij}\in\mathcal{T}_{i}$,
$t_{j^{\prime}}\in\mathcal{T}_{m}$ and $\delta_{jj^{\prime}}=0$, otherwise.

The proofs of Theorem 1 and the following Corollary are included in Appendix
A.

###### Corollary 1

The $\boldsymbol{\mu}_{k}$, $\boldsymbol{\nu}_{k}$, $\boldsymbol{\alpha}_{i}$
and $\boldsymbol{\gamma}_{i}$ can be obtained through a backfitting algorithm
or the Gauss-Seidel method, iterating the following two steps until
convergence:

  1. (a) for each $i$, $\boldsymbol{\hat{\alpha}}_{i}=(\boldsymbol{\phi}_{\alpha_{i}}^{\prime}\boldsymbol{S}_{U_{i}}^{-1}\boldsymbol{\phi}_{\alpha_{i}})^{-1}\boldsymbol{\phi}_{\alpha_{i}}^{\prime}\boldsymbol{S}_{U_{i}}^{-1}\boldsymbol{\tilde{Y}}_{i}$ and $\boldsymbol{\hat{\gamma}}_{i}=\boldsymbol{S}_{U_{i}}^{-1}\left\\{\boldsymbol{I}-\boldsymbol{\phi}_{\alpha_{i}}(\boldsymbol{\phi}_{\alpha_{i}}^{\prime}\boldsymbol{S}_{U_{i}}^{-1}\boldsymbol{\phi}_{\alpha_{i}})^{-1}\boldsymbol{\phi}_{\alpha_{i}}^{\prime}\boldsymbol{S}_{U_{i}}^{-1}\right\\}\boldsymbol{\tilde{Y}}_{i}$, where $\boldsymbol{S}_{U_{i}}=\boldsymbol{R}_{U_{i}1}+n_{i}\lambda_{U_{i}}\boldsymbol{I}$ and $\boldsymbol{\tilde{Y}}_{i}=\boldsymbol{Y}_{i}-\boldsymbol{\Delta}_{i}\boldsymbol{\phi}_{\mu}\boldsymbol{\hat{\mu}}_{k_{i}}-\boldsymbol{\Delta}_{i}\boldsymbol{R}_{M_{1}}\boldsymbol{\hat{\nu}}_{k_{i}}$;

  2. (b) for each $k$, $\boldsymbol{\hat{\mu}}_{k}=(\boldsymbol{\phi}_{\mu}^{\prime}\boldsymbol{\Delta}^{\prime}\boldsymbol{S}_{M_{k}}^{-1}\boldsymbol{\Delta}\boldsymbol{\phi}_{\mu})^{-1}\boldsymbol{\phi}_{\mu}^{\prime}\boldsymbol{\Delta}^{\prime}\boldsymbol{S}_{M_{k}}^{-1}\boldsymbol{\tilde{Y}}_{k}$ and
$\boldsymbol{\hat{\nu}}_{k}=\boldsymbol{S}_{M_{k}}^{-1}\left\\{\boldsymbol{I}-\boldsymbol{\Delta}\boldsymbol{\phi}_{\mu}(\boldsymbol{\phi}_{\mu}^{\prime}\boldsymbol{\Delta}^{\prime}\boldsymbol{S}_{M_{k}}^{-1}\boldsymbol{\Delta}\boldsymbol{\phi}_{\mu})^{-1}\boldsymbol{\phi}_{\mu}^{\prime}\boldsymbol{\Delta}^{\prime}\boldsymbol{S}_{M_{k}}^{-1}\right\\}\boldsymbol{\tilde{Y}}_{k}$,
where
$\boldsymbol{S}_{M_{k}}=\boldsymbol{\Delta}\boldsymbol{R}_{M_{1}}+\lambda_{M_{k}}\boldsymbol{I}$,
$\boldsymbol{\tilde{Y}}_{k}=\displaystyle\sum_{i:k_{i}=k}\frac{1}{n_{i}}\boldsymbol{\Delta}_{i}^{\prime}\left(\boldsymbol{Y}_{i}-\boldsymbol{\phi}_{\alpha_{i}}\boldsymbol{\hat{\alpha}}_{i}-\boldsymbol{R}_{U_{i1}}\boldsymbol{\hat{\gamma}}_{i}\right)$
and
$\boldsymbol{\Delta}=\displaystyle\sum_{i:k_{i}=k}\frac{1}{n_{i}}\boldsymbol{\Delta}_{i}^{\prime}\boldsymbol{\Delta}_{i}$.

##  3 Posterior Computation

Although we can obtain $\hat{M}_{k}(t)$ and $\hat{U}_{i}(t)$ by the
backfitting algorithm outlined in Corollary 1 and estimate $\lambda_{M_{k}}$
and $\lambda_{U_{i}}$ through generalized cross validation (Chap. 4, Wahba,
1990), it is unclear how to allow $\lambda_{U_{i}}$ to depend on covariates.
In addition, when $n$ is large, it is computational infeasible to invert the
$n\times n$ matrix $\boldsymbol{S}_{M_{k}}$ involved in the backfitting
algorithm. Instead, we propose a Markov chain Monte Carlo (MCMC) algorithm for
posterior computation that solves these problems. The algorithm achieves
computational efficiency by leveraging on the Markovian property of SDEs and
samples $M_{k}(t)$ and $U_{i}(t)$ through the simulation smoother (Durbin and
Koopman, 2002), which requires the following Proposition.

###### Proposition 1

Let $X(t)$ denote a (r-1)th-order integral Wiener process, defined by the
stochastic differential equation $D^{r}X(t)=\dot{W}(t)$. Consequently, the
$\boldsymbol{X}_{j}=\\{X(t_{j}),D^{1}X(t_{j}),\cdots,D^{r-1}X(t_{j})\\}^{\prime}$,
$j=1,2,\cdots,n$, follows a state equation

$\boldsymbol{X}_{j+1}=\boldsymbol{G}_{j}\boldsymbol{X}_{j}+\boldsymbol{\omega}_{j},$

where
$\boldsymbol{G}_{j}=\sum_{k=0}^{r}\frac{\delta_{j}^{k}}{k!}\boldsymbol{C}^{k}$
and
$\boldsymbol{\omega}_{j}\sim\mathsf{N}_{r}(\boldsymbol{0},\boldsymbol{W}_{j})$
with
$\boldsymbol{W}_{j}=\int_{0}^{\delta_{j}}\exp\\{\boldsymbol{C}(\delta_{j}-u)\\}\boldsymbol{D}\boldsymbol{D}^{\prime}\exp\\{\boldsymbol{C}^{\prime}(\delta_{j}-u)\\}du$,
$\boldsymbol{C}=(c_{ll^{\prime}})_{r\times r}$, $c_{ll^{\prime}}=1$ when
$l^{\prime}=l+1$ and $c_{ll^{\prime}}=0$ otherwise,
$\boldsymbol{D}=(0,0,\cdots,1)^{\prime}$ and $\delta_{j}=t_{j+1}-t_{j}$.

The proof is in Appendix A. Finally, we outline the proposed MCMC as follows.

  1. (1) Given $M_{k_{i}}(t_{ij})$, $\sigma^{2}_{\varepsilon}$ and $\sigma^{2}_{U_{i}}$, sample $U_{i}(t_{ij})$, $i=1,2,\cdots,m$, $j=0,1,\cdots,n_{i}$. Let $Y_{U_{ij}}=Y_{i}(t_{ij})-M_{k_{i}}(t_{ij})$ and the SVR model for the $i$th subject can be expressed as the following state space model (Jones, 1993; Durbin and Koopman, 2001), from which we can draw samples of $U_{i}(t_{ij})$ and its derivatives using the simulation smoother.

$\displaystyle Y_{U_{ij}}$
$\displaystyle=\boldsymbol{F}_{U_{ij}}\boldsymbol{U}_{{ij}}+\varepsilon_{U_{ij}},$
$\displaystyle\boldsymbol{U}_{{i(j+1)}}$
$\displaystyle=\boldsymbol{G}_{U_{ij}}\boldsymbol{U}_{{ij}}+\sigma_{U_{i}}\boldsymbol{\omega}_{U_{ij}},$

where $\boldsymbol{F}_{U_{ij}}=(1,0,\cdots,0)$,
$\boldsymbol{U}_{{ij}}=\\{U_{i}(t_{ij}),D^{1}U_{i}(t_{ij}),\cdots,D^{q-1}U_{i}(t_{ij})\\}^{\prime}$
and
$\varepsilon_{U_{ij}}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\mathsf{N}_{1}(0,\sigma^{2}_{\varepsilon})$.
Similar to the $\boldsymbol{G}_{j}$, $\boldsymbol{\omega}_{j}$ and
$\boldsymbol{W}_{j}$ in Proposition 1, the $\boldsymbol{G}_{U_{ij}}$,
$\boldsymbol{\omega}_{U_{ij}}$ and $\boldsymbol{W}_{U_{ij}}$ follow the same
specifications with $r=q$.

  2. (2) Given $U_{i}(t_{j})$, $\sigma^{2}_{\varepsilon}$ and $\sigma^{2}_{M_{k}}$, sample $M_{k}(t_{j})$, $k=1,2,\cdots,g$, $j=0,1,\cdots,n$. Similarly, we rewrite the SVR model for the $k$th group as the following state space model and then sample $M_{k_{i}}(t_{ij})$ and its derivatives by the simulation smoother.

$\displaystyle\boldsymbol{Y}_{M_{kj}}$
$\displaystyle=\boldsymbol{F}_{M_{kj}}\boldsymbol{M}_{kj}+\boldsymbol{\varepsilon}_{M_{kj}},$
$\displaystyle\boldsymbol{M}_{{k(j+1)}}$
$\displaystyle=\boldsymbol{G}_{M_{kj}}\boldsymbol{M}_{{kj}}+{\sigma}_{M_{k}}\boldsymbol{\omega}_{M_{kj}},$

where $\boldsymbol{Y}_{M_{kj}}=(Y_{M_{kj}}^{i})_{m\times 1}$,
$\boldsymbol{M}_{{kj}}=\\{M_{k}(t_{j}),D^{1}M_{k}(t_{j}),\cdots,D^{p-1}M_{k}(t_{j})\\}^{\prime}$,
$\boldsymbol{F}_{M_{kj}}=(F_{M_{kj}}^{il})_{m\times p}$ and
$\boldsymbol{\varepsilon}_{M_{kj}}=\text{diag}(\varepsilon_{M_{kj}}^{1},\varepsilon_{M_{kj}}^{2},\cdots,\varepsilon_{M_{kj}}^{m})$.
When $i$th subject has an observation at time $t_{j}$ and $k_{i}=k$,
$Y_{M_{kj}}^{i}=Y_{i}(t_{j})-U_{i}(t_{j})$, $F_{M_{kj}}^{i1}=1$ and
$\varepsilon_{M_{kj}}^{i}\sim\mathsf{N}_{1}(0,\sigma^{2}_{\varepsilon})$.
Otherwise, $Y_{M_{kj}}^{i}=F_{M_{kj}}^{il}=\varepsilon_{M_{kj}}^{i}=0$. The
$\boldsymbol{G}_{M_{kj}}$, $\boldsymbol{\omega}_{M_{kj}}$ and
$\boldsymbol{W}_{M_{kj}}$ are given by Proposition 1 with $r=p$.

  3. (3a) Given $M_{k_{i}}(t_{ij})$ and $U_{i}(t_{ij})$, $i=1,2,\cdots,m$, $j=1,2,\cdots,n_{i}$, sample $\sigma^{2}_{\varepsilon}$ from the posterior distribution $\textsf{invGamma}\left(a+\frac{1}{2}\sum_{i=1}^{m}n_{i},b+\frac{1}{2}\sum_{i=1}^{m}\sum_{j=1}^{n_{i}}\left\\{Y_{i}(t_{ij})-M_{k_{i}}(t_{ij})-U_{i}(t_{ij})\right\\}^{2}\right)$.

  4. (3b) Given $\boldsymbol{U}_{{i0}}$, sample $\sigma^{2}_{U_{0}}$ from the posterior distribution $\textsf{invGamma}\left(a+\frac{mq}{2},b+\frac{1}{2}\sum_{i=0}^{m}\boldsymbol{U}_{{i0}}^{\prime}\boldsymbol{U}_{{i0}}\right)$.

  5. (3c) Given $\boldsymbol{M}_{{kj}}$, sample $\sigma_{M_{k}}^{2}$ from the posterior distribution
$\textsf{invGamma}\left(a+\frac{np}{2},b+\frac{1}{2}\sum_{j=0}^{n-1}(\boldsymbol{M}_{{k(j+1)}}-\boldsymbol{G}_{M_{kj}}\boldsymbol{M}_{{kj}})^{\prime}\boldsymbol{W}_{M_{kj}}^{-1}(\boldsymbol{M}_{{k(j+1)}}-\boldsymbol{G}_{M_{kj}}\boldsymbol{M}_{{kj}})\right)$.

  6. (3d) Given $\boldsymbol{U}_{{ij}}$, $\boldsymbol{\beta}$ and $\sigma^{2}$, sample $\sigma^{2}_{U_{i}}$ using a Metropolis-Hasting algorithm. We choose $\sigma^{2}_{U_{i}}\sim\textsf{invGamma}(a,b)$ as the proposal prior distribution and a proposal $\sigma^{2*}_{U_{i}}$ can be easily drawn from $\textsf{invGamma}\left(a+\frac{n_{i}q}{2},b+\frac{1}{2}\sum_{j=0}^{n_{i}-1}(\boldsymbol{U}_{{i(j+1)}}-\boldsymbol{G}_{U_{ij}}\boldsymbol{U}_{{ij}})^{\prime}\boldsymbol{W}_{U_{ij}}^{-1}(\boldsymbol{U}_{{i(j+1)}}-\boldsymbol{G}_{U_{ij}}\boldsymbol{U}_{{ij}})\right)$ the corresponding proposal posterior distribution. The $\sigma^{2*}_{U_{i}}$ will be accepted with the following probability and discarded otherwise with $\sigma^{2}_{U_{i}}$ unchanged,

$\min\left\\{\frac{f_{\mathsf{LN}}(\sigma^{2*}_{U_{i}}\mid\boldsymbol{x}^{\prime}_{i}\boldsymbol{\beta},\sigma^{2})\prod_{j=0}^{n_{i}-1}f_{\mathsf{N_{q}}}(\boldsymbol{U}_{{i(j+1)}}-\boldsymbol{G}_{U_{ij}}\boldsymbol{U}_{{ij}}\mid\boldsymbol{0},\sigma^{2*}_{U_{i}}\boldsymbol{W}_{U_{ij}})f_{\mathsf{iG}}(\sigma^{2}_{U_{i}}\mid
a_{U_{i}},b_{U_{i}})}{f_{\mathsf{LN}}(\sigma^{2}_{U_{i}}\mid\boldsymbol{x}^{\prime}_{i}\boldsymbol{\beta},\sigma^{2})\prod_{j=0}^{n_{i}-1}f_{\mathsf{N_{q}}}(\boldsymbol{U}_{{i(j+1)}}-\boldsymbol{G}_{U_{ij}}\boldsymbol{U}_{{ij}}\mid\boldsymbol{0},\sigma^{2}_{U_{i}}\boldsymbol{W}_{U_{ij}})f_{\mathsf{iG}}(\sigma^{2*}_{U_{i}}\mid
a_{U_{i}},b_{U_{i}})},1\right\\},$

where $f_{\mathsf{LN}}$, $f_{\mathsf{N_{q}}}$ and $f_{\mathsf{iG}}$ denote the
log-normal, $q$-dimensional normal and inverse gamma probability density
functions respectively with $a_{U_{i}}=a+\frac{n_{i}q}{2}$,
$b_{U_{i}}=b+\frac{1}{2}\sum_{j=0}^{n_{i}-1}(\boldsymbol{U}_{{i(j+1)}}-\boldsymbol{G}_{U_{ij}}\boldsymbol{U}_{{ij}})^{\prime}\boldsymbol{W}_{U_{ij}}^{-1}(\boldsymbol{U}_{{i(j+1)}}-\boldsymbol{G}_{U_{ij}}\boldsymbol{U}_{{ij}})$.

  7. (4) Given $\sigma^{2}_{U_{i}}$, sample $\boldsymbol{\beta}$ and $\sigma^{2}$. Let $\boldsymbol{Z}=(\log\sigma^{2}_{U_{1}},\log\sigma^{2}_{U_{2}},\cdots,\log\sigma^{2}_{U_{m}})^{\prime}$, $\boldsymbol{\hat{\beta}}=(\boldsymbol{X}^{\prime}\boldsymbol{X})^{-1}\boldsymbol{X}^{\prime}\boldsymbol{Z}$ and $\hat{\sigma}^{2}=\frac{(\boldsymbol{Z}-\boldsymbol{X}\boldsymbol{\hat{\beta}})^{\prime}(\boldsymbol{Z}-\boldsymbol{X}\boldsymbol{\hat{\beta}})}{m-k}$. We draw $\tau$ from Chi-squared distribution with $m-k$ degrees of freedom and set $\sigma^{2}=\frac{(m-k)\hat{\sigma}^{2}}{\tau}$ and then sample $\boldsymbol{\beta}$ from $\textsf{N}_{m}\left(\boldsymbol{\hat{\beta}},\sigma^{2}(\boldsymbol{X}^{\prime}\boldsymbol{X})^{-1}\right)$.

##  4 Simulation

We carry out two simulation studies to evaluate the performance of the
proposed method and compare it to alternative methods including nature cubic
spline (NCS, Wahba, 1990) and functional principal components analysis (FPCA,
Yao et al., 2005). The comparison focuses on performance in estimating the
trajectory $M_{k_{i}}(t)+U_{i}(t)$, the volatility $\sigma^{2}_{U_{i}}$ and
the coefficients $\boldsymbol{\beta}$.

The first simulation study is designed to investigate the consequence of
ignoring either similarity or heterogeneity of volatilities when they are
present. One hundred replicated datasets, each consisting of 100 trajectories,
are sampled from the SVR model, in which the log-transformed volatilities are
varying and normally distributed. More precisely, we choose
$\boldsymbol{\beta}=(0,0.6,2)^{\prime}$ and
$\boldsymbol{x}_{i}=(1,x_{i1},x_{i2})^{\prime}$ with $x_{i1}$ and $x_{i2}$
sampled from
$x_{i1}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\mathsf{Bernoulli}(0.4)$
and
$x_{i2}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\mathsf{N}_{1}(0,0.25)$
respectively. Given $\boldsymbol{\beta}$ and $\boldsymbol{x}_{i}$,
volatilities $\sigma^{2}_{U_{i}}$’s can be drawn from
$\log(\sigma^{2}_{U_{i}})\sim\mathsf{N}_{1}(\boldsymbol{x}^{\prime}_{i}\boldsymbol{\beta},1)$.
Along with $\sigma^{2}_{M_{1}}=\sigma^{2}_{M_{2}}=10$,
$\sigma^{2}_{\varepsilon}=1$, $p=2$ and $q=1$, $M_{1}(t)$, $M_{2}(t)$,
$U_{i}(t)$ and $\varepsilon_{i}(t)$ are sampled at
$t\in\\{0.2,0.4,\cdots,4\\}$ from equations (2) and (3) and the distribution
of measurement error $\varepsilon_{i}(t)$. Twenty percent of samples are
removed completely at random, resulting in on average 16 unequally spaced
observations per subject. Finally, the $i$th subject is randomly assigned to
one of the two groups with equal probability and $Y_{i}(t)$ is obtained from
observation equation (1).

The simulated datasets are analyzed by three methods, SVR, NCS and FPCA. We
first apply the SVR approach, using the proposed MCMC algorithm with 15,000
iterations and keeping every 5th of the last 10,000 samples for posterior
analysis. It takes about 80 minutes on a PC with 2.33GHz Intel(R) Xeon(R) CPU.
Posterior means are chosen as the estimates of $M_{k_{i}}(t)+U_{i}(t)$,
$\sigma^{2}_{U_{i}}$ and $\boldsymbol{\beta}$. Additionally, the trajectories
$M_{k_{i}}(t)+U_{i}(t)$’s are estimated by NCS for one subject at a time, and
by FPCA for all subjects within a group and separately by the group, taking
about 1 minute and 2 minutes on the same PC respectively. For NCS and FPCA
methods, we may also estimate covariate effects on volatility through a two-
stage method: estimating empirical volatility by
$\frac{1}{n_{i}}\sum_{j=1}^{n_{i}-1}\frac{(\hat{U}_{i,j+1}-\hat{U}_{i,j})^{2}}{t_{i,j+1}-t_{i,j}}$
in the first stage with $\hat{U}_{i,j}$ the estimate of $U_{i}(t)$ at time
$t_{ij}$; and in the second stage, empirical volatilities are regressed on
covariates to obtain the estimate of $\boldsymbol{\beta}$.

For each simulated dataset, we calculate average squared error (ASE) for the
trajectory
$\text{ASE}(M+U)=\frac{1}{m}\sum_{i=1}^{m}\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\left\\{\hat{M}_{k_{i}}(t_{ij})+\hat{U}_{i}(t_{ij})-M_{k_{i}}(t_{ij})-U_{i}(t_{ij})\right\\}^{2}$,
ASE for log volatility
$\text{ASE}\\{log(\sigma^{2}_{U})\\}=\frac{1}{m}\sum_{i=1}^{m}\left\\{log(\hat{\sigma}^{2}_{U_{i}})-log(\sigma^{2}_{U_{i}})\right\\}^{2}$,
and squared errors (SE) for coefficient estimates
$\text{SE}(\beta_{l})=(\hat{\beta}_{l}-\beta_{l})^{2}$, $l=0,1,2$. Table 1
reports means of $\text{ASE}(M+U)$, $\text{ASE}\\{log(\sigma^{2}_{U})\\}$ and
$\text{SE}(\beta_{l})$ across 100 replicate datasets. MASEs and MSEs by NCS
and FPCA approaches are significantly inflated, for example, being doubled and
tripled in $\text{MASE}(M+U)$ respectively, compared to SVR. We randomly
select a data set and take a close examination. We calculate the individual
ASE of the trajectory
$\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\left\\{\hat{M}_{k_{i}}(t_{ij})+\hat{U}_{i}(t_{ij})-M_{k_{i}}(t_{ij})-U_{i}(t_{ij})\right\\}^{2}$
and select the top four subjects with the largest individual ASEs with respect
to NCS and FPCA approaches respectively.

Figure 2 shows estimates of the trajectory for six subjects. The plots reveal
that the increased MASEs or MSEs by NCS and FPCA are caused by different
reasons. NCS approach, treating one trajectory a time, ignores the similarity
between the subjects within a group, leading to over fitting true trajectories
(e.g. Figure 2(d) and 2(e)) with both over and under estimated volatilities.
On the other hand, FPCA approach omits the heterogeneity of the subjects
within a group; inflated $\text{MASE}(M+U)$ are mainly contributed by a few
subjects whose trajectory fluctuates with significantly higher volatility but
is overly smoothed (e.g. Figure 2(b) and 2(d)); and under the assumption of
similar smoothness, the estimates of volatility are largely under estimated.
Although this simulation study is in favor of SVR approach by design, the
scenario we consider is nevertheless realistic in practice and the simulation
results reveal the drawbacks of omitting similarity or heterogeneity of
volatilities by alternative approaches.

Our second simulation study is conducted to evaluate the performance of SVR,
NCS and FPCA when volatilities are homogeneous with no covariate effects. As
in the first simulation study, 100 replicate datasets are generated, each
consisting of 100 trajectories at $t\in\\{0.2,0.4,\cdots,4\\}$; twenty percent
of data points are deleted completely at random; and subjects are assigned to
one of two groups with equal probability. The observations are generated from
$Y_{i}(t)=10\\{t+sin(t)\\}+0.6\alpha_{1i}cos(\pi t/10)+0.2\alpha_{2i}sin(\pi
t/10)+\varepsilon_{i}(t)$ for subjects in the first group and from
$Y_{i}(t)=10\\{t+cos(t)\\}+0.5\alpha_{1i}cos(\pi t/10)+0.3\alpha_{2i}sin(\pi
t/10)+\varepsilon_{i}(t)$ for the ones in the second group, with
$\alpha_{1i}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\mathsf{N}_{1}(0,4)$,
$\alpha_{2i}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\mathsf{N}_{1}(0,1)$
and
$\varepsilon_{i}(t)\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\mathsf{N}_{1}(0,1)$.
The SVR, NCS and FPCA approaches are applied and $\text{MASE}(M+U)$ is
presented in Table 1, in which SVR and FPCA approaches show close
$\text{MASE}(M+U)$, both smaller than the one by NCS approach. This suggests
that SVR is no worse than FPCA for the cases with homogeneous volatilities.

In short, through the two simulation studies, we demonstrate that SVR achieves
substantial gains in terms of reducing the ASEs or SEs of the estimates of the
trajectory, volatility and covariate effect when volatilities are
heterogeneous, and works at least as well as FPCA approach when volatilities
are homogeneous.

Table 1: The mean of squared errors or average square errors of the estimates of trajectory, volatility and covariate effect across 100 replicate datasets . |  |  | Case I |  |  | Case II
---|---|---|---|---|---|---
method | $M+U$ | $\log(\sigma^{2}_{U})$ | $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$ | $M+U$
SVR | 0.345 | 0.614 | 0.043 | 0.081 | 0.075 | 1.122
NCS | 0.609 | 1.297 | 0.089 | 0.165 | 1.724 | 1.477
FPCA | 1.099 | 2.966 | 1.144 | 0.185 | 1.969 | 1.185

(a) ID=17

(b) ID=30

(c) ID=31

(d) ID=39

(e) ID=44

(f) ID=86

Figure 2:  The plots of observation ($\circ$) and trajectory at time $t_{ij}$
($\times$), as well as estimates of trajectory $M_{k_{i}}(t)+U_{i}(t)$ by SVR
(—), NCS ($---$) and FPCA ($\cdotp\cdotp\cdotp$) approaches, for six subjects
in one simulated dataset with the largest individual average squared errors
$\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\left\\{\hat{M}_{k_{i}}(t_{ij})+\hat{U}_{i}(t_{ij})-M_{k_{i}}(t_{ij})-U_{i}(t_{ij})\right\\}^{2}$.

##  5 Applications

It is a common practice to monitor the blood pressure of pregnant woman during
pregnancy. Despite the trend of blood pressure being well studied, its
fluctuation has been little addressed and the factors associated with the
fluctuation are largely unknown. In this analysis, we apply the proposed SVR
approach to analyze longitudinal blood pressure measurements in HPHB study,
aiming to investigate the stability of blood pressure trajectories and
identify the associated factors. The data set consists of 106 non-Hispanic
white (NHW) and 176 non-Hispanic black (NHB) women whose first blood pressure
measurement is collected before the 16th week of gestation and the last one no
earlier than the 37th week of gestation. Most of subjects have 9 ($35.10\%$ of
them), 10 ($29.28\%$) or 11 ($14.98\%$) measurements spaced at irregularly
times. The covariates we focused on include race as NHW vs NHB, mother’s age
group as age $>$ 35 vs age $\leq$ 35, obesity as yes vs no, preeclampsia as
yes vs no, parity group as parity $>$ 0 vs parity $=$ 0, and smoking as yes vs
no. We run the proposed MCMC algorithm for 15,000 iterations, discard the
first 5,000 and retain every 5th of the remaining samples for posterior
analysis. The trace plots and autocorrelation plots suggest the algorithm
converges fast and mixes well. Posterior summary of selected parameters is
presented in Table 2.

The panels LABEL:sub@fig:MAPMNHW and LABEL:sub@fig:MAPMNHB of Figure 3 show
posterior means and $95\%$ credible intervals of the average blood pressure
for NHW and NHB groups respectively, which share a common pattern: decreasing
till the late stage of the second trimester (about 20 to 25 weeks) and then
increasing toward the pre-pregnancy level. Within the ethnic group,
significant heterogeneity exists in terms of the stability of the blood
pressure trajectory at the individual level. As Figure 3
LABEL:sub@fig:MAPVolPost indicated, posterior means of volatility vary from
-0.5 to 2 in the logarithmic scale, suggesting some subjects’ trajectories are
parallel to the group blood pressure trajectory with very small volatilities
while others may significantly depart from the group blood pressure
trajectory.

Most interesting, we find that obesity and preeclampsia are associated with
blood pressure volatility, with their $95\%$ credible intervals not covering
zero in Figure 3 LABEL:sub@fig:MAPCoeff. This implies that pregnant women with
obesity and/or preeclampsia are more likely to demonstrate irregular patterns
of blood pressure relative to their ethnic group. We further examine the
characteristics of these subjects with extreme volatilities (results not
shown). Among the top eight subjects presenting with the largest volatilities,
most of them are NHB with obesity and preeclampsia, do not smoke and give
birth to a baby for the first time; half of them are younger than 35. For the
eight subjects with the smallest volatilities, they are surprising
homogeneous, i.e. all of them being NHW (except one) without obesity and
preeclampsia, younger than 35, not smoking and giving birth to a baby before.

Table 2: Blood pressure data: Posterior summary of parameters in the SVR model. Parameter | Mean | Mode | SD | 95% HPD invteral
---|---|---|---|---
$\sigma^{2}_{\varepsilon}$ | 17.807 | 17.818 | 0.694 | [16.389, 19.106]
$\sigma^{2}_{M_{1}}$ | 0.236 | 0.187 | 0.181 | [0.040, 0.556]
$\sigma^{2}_{M_{2}}$ | 0.204 | 0.162 | 0.148 | [0.042, 0.472]
$\sigma^{2}_{U_{0}}$ | 46.729 | 46.412 | 4.776 | [38.263, 56.619]
$\sigma^{2}$ | 0.734 | 0.741 | 0.333 | [0.082, 1.295]

(a)

(b)

(c)

(d)

Figure 3:  The posterior means and $95\%$ highest posterior density (HPD)
credible intervals for (a) the blood pressure during the 2nd and 3rd
trimesters for non-Hispanic white group; (b) the blood pressure during the 2nd
and 3rd trimesters for non-Hispanic black group; (c) the volatility in the
logarithmic scale; (d) covariate effects.

##  6 Discussion

We have proposed a stochastic volatility regression model to investigate the
volatility and its association with covariates for multi-subject functional
data. As an important dynamic feature, volatility measures the stability of
the biological process. The analysis of volatility not only reveals its
heterogeneity among subjects but also its dependence on the covariates of
interest. Complementing the current FDA methods which mainly focus on the
trend of trajectory and its derivatives, the SVR method initiates the
exploration of stability of functional data. As illustrated with the blood
pressure data, our view is that substantial new insights can be obtained in a
rich variety of biomedical applications by studying volatility.

## Acknowledgments

This work was supported by Award Number R01ES017436 and R01ES17240 from the
National Institute of Environmental Health Sciences, by funding from the
National Institutes of Health (5P2O-RR020782-O3) and the U.S. Environmental
Protection Agency (RD-83329301-0), and by the Intramural Research Program of
the National Cancer Institute, National Institutes of Health. The content is
solely the responsibility of the authors and does not necessarily represent
the official views of the National Institute of Environmental Health Sciences,
the National Institutes of Health or the U.S. Environmental Protection Agency.

##  Appendix A Appendix: Proofs of Theoretical Results

###  A.1 Proof of Theorem 1

By the RKHS theory of the polynomial smoothing spline (Wahba, 1990, Section
1.2), we have

$\displaystyle{M}_{k}(t)$
$\displaystyle={M}_{k0}(t)+{M}_{k1}(t)=\sum_{l=0}^{p-1}\mu_{kl}\phi_{l}(t)+\sum_{j=1}^{n}\nu_{kj}\mathcal{R}_{M_{1}}(t_{j},t)+\eta_{M_{1}}(t),$

and

$\int_{\mathcal{T}}\left\\{D^{p}M_{k}(t)\right\\}^{2}dt=\boldsymbol{\nu}^{\prime}_{k}\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k}+\langle\eta_{M_{1}}(\cdot),\eta_{M_{1}}(\cdot)\rangle_{\mathcal{H}_{\mathcal{R}_{M_{1}}}},$

where $M_{0}(t)\in\mathcal{H}_{\mathcal{R}_{M_{0}}}$ and
$M_{1}(t)\in\mathcal{H}_{\mathcal{R}_{M_{1}}}$, the RKHSs
$\mathcal{H}_{\mathcal{R}_{M_{0}}}=\left\\{f(t):D^{p}f(t)=0,t\in\mathcal{T}\right\\}$
and $\mathcal{H}_{\mathcal{R}_{M_{1}}}=\left\\{f(t):D^{i}f(t)\text{ absolutely
continuous for
}i=0,1,\cdots,p-1,D^{p}f(t)\in\mathcal{L}_{2}(\mathcal{T})\right\\}$;
$\mathcal{R}_{M_{0}}(s,t)$ and $\mathcal{R}_{M_{1}}(s,t)$ are reproducing
kernels defined in Lemma 1;
$\eta_{M_{1}}(\cdot)\in\mathcal{H}_{\mathcal{R}_{M_{1}}}$ is orthogonal to
$\mathcal{R}_{M_{1}}(t_{j},\cdot)$ with inner product
$\langle\mathcal{R}_{M_{1}}(t_{j},\cdot),\eta_{M_{1}}(\cdot)\rangle_{\mathcal{H}_{\mathcal{R}_{M_{1}}}}=\int_{\mathcal{T}}D^{p}\mathcal{R}_{M_{1}}(t_{j},u)D^{p}\eta_{M_{1}}(u)du=0$
for $j=1,2,\cdots,J$;
$\mathcal{L}_{2}(\mathcal{T})=\left\\{f(t):\int_{\mathcal{T}}f^{2}(t)dt<\infty\right\\}$
is the space of squared integrable functions defined on index set
$\mathcal{T}$.

Similarly,

$\displaystyle{U}_{i}(t)$
$\displaystyle={U}_{i0}(t)+{U}_{i1}(t)=\sum_{l=0}^{q-1}\alpha_{il}\phi_{l}(t)+\sum_{j=1}^{n_{i}}\gamma_{ij}\mathcal{R}_{U_{1}}(t_{ij},t)+\eta_{U_{i1}}(t),$

and

$\int_{\mathcal{T}}\left\\{D^{q}U_{i}(t)\right\\}^{2}dt=\boldsymbol{\gamma}_{i}^{\prime}\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i}+\langle\eta_{U_{i1}}(\cdot),\eta_{U_{i1}}(\cdot)\rangle_{\mathcal{H}_{\mathcal{R}_{U_{1}}}},$

where $U_{i0}(t)\in\mathcal{H}_{\mathcal{R}_{U_{0}}}$ and
$U_{i1}(t)\in\mathcal{H}_{\mathcal{R}_{U_{1}}}$ with
$\mathcal{H}_{\mathcal{R}_{U_{0}}}=\left\\{f(t):D^{q}f(t)=0,t\in\mathcal{T}\right\\}$
and $\mathcal{H}_{\mathcal{R}_{U_{1}}}=\left\\{f(t):D^{i}f(t)\text{ absolutely
continuous for
}i=0,1,\cdots,q-1,D^{q}f(t)\in\mathcal{L}_{2}(\mathcal{T})\right\\}$ the RKHSs
with reproducing kernel $\mathcal{R}_{U_{0}}(s,t)$ and
$\mathcal{R}_{U_{1}}(s,t)$ defined in Lemma 1;
$\eta_{U_{i1}}(\cdot)\in\mathcal{H}_{\mathcal{R}_{U_{1}}}$ is orthogonal to
$\mathcal{R}_{U_{1}}(t_{ij},\cdot)$ for $j=1,2,\cdots,n_{i}$.

Hence, the double-penalized sum-of-squares (4) can be written as

$\displaystyle\textsf{DPSS}=$
$\displaystyle\sum_{i=1}^{m}\frac{1}{n_{i}}(\boldsymbol{Y}_{i}-\boldsymbol{\Delta}_{i}\boldsymbol{\phi}_{\mu}\boldsymbol{\mu}_{k_{i}}-\boldsymbol{\Delta}_{i}\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k_{i}}-\boldsymbol{\phi}_{\alpha_{i}}\boldsymbol{\alpha}_{i}-\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i})^{\prime}\times$
$\displaystyle(\boldsymbol{Y}_{i}-\boldsymbol{\Delta}_{i}\boldsymbol{\phi}_{\mu}\boldsymbol{\mu}_{k_{i}}-\boldsymbol{\Delta}_{i}\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k_{i}}-\boldsymbol{\phi}_{\alpha_{i}}\boldsymbol{\alpha}_{i}-\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i})+$
$\displaystyle\sum_{k=1}^{g}\lambda_{M_{k}}\boldsymbol{\nu}^{\prime}_{k}\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k}+\sum_{i=1}^{m}\lambda_{U_{i}}\boldsymbol{\gamma}^{\prime}_{i}\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i}+$
$\displaystyle\sum_{k=1}^{g}\lambda_{M_{k}}\langle\eta_{M_{1}}(\cdot),\eta_{M_{1}}(\cdot)\rangle_{\mathcal{H}_{\mathcal{R}_{M_{1}}}}+\sum_{i=1}^{m}\lambda_{U_{i}}\langle\eta_{U_{i1}}(\cdot),\eta_{U_{i1}}(\cdot)\rangle_{\mathcal{H}_{\mathcal{R}_{U_{1}}}},$

which is minimized at
$\langle\eta_{M_{1}}(\cdot),\eta_{M_{1}}(\cdot)\rangle_{\mathcal{H}_{\mathcal{R}_{M_{1}}}}=\langle\eta_{U_{i1}}(\cdot),\eta_{U_{i1}}(\cdot)\rangle_{\mathcal{H}_{\mathcal{R}_{U_{1}}}}=0$,
leading to $\eta_{M_{1}}(\cdot)=\eta_{U_{i1}}=0$.

###  A.2 Proof of Corollary 1

Taking partial derivatives of double penalized sum-of-squares in Corollary 1
with respective to $\boldsymbol{\mu}_{k}$, $\boldsymbol{\nu}_{k}$,
$\boldsymbol{\alpha}_{i}$ and $\boldsymbol{\gamma}_{i}$ respectively and
setting them to zeros, we have

$\displaystyle\frac{\partial\;\textsf{DPSS}}{\partial\;\boldsymbol{\mu}_{k}}$
$\displaystyle=\displaystyle\sum_{i:k_{i}=k}\frac{1}{n_{i}}\boldsymbol{\phi}_{\mu}^{\prime}\boldsymbol{\Delta}_{i}^{\prime}(\boldsymbol{\Delta}_{i}\boldsymbol{\phi}_{\mu}\boldsymbol{\mu}_{k_{i}}+\boldsymbol{\Delta}_{i}\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k_{i}}+\boldsymbol{\phi}_{\alpha_{i}}\boldsymbol{\alpha}_{i}+\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i}-\boldsymbol{Y}_{i})=\boldsymbol{0},$
$\displaystyle\frac{\partial\;\textsf{DPSS}}{\partial\;\boldsymbol{\nu}_{k}}$
$\displaystyle=\displaystyle\sum_{i:k_{i}=k}\frac{1}{n_{i}}\boldsymbol{R}_{M_{1}}\boldsymbol{\Delta}_{i}^{\prime}(\boldsymbol{\Delta}_{i}\boldsymbol{\phi}_{\mu}\boldsymbol{\mu}_{k_{i}}+\boldsymbol{\Delta}_{i}\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k_{i}}+\boldsymbol{\phi}_{\alpha_{i}}\boldsymbol{\alpha}_{i}+\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i}-\boldsymbol{Y}_{i})+\lambda_{M_{k}}\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k}=\boldsymbol{0},$
$\displaystyle\frac{\partial\;\textsf{DPSS}}{\partial\;\boldsymbol{\alpha}_{i}}$
$\displaystyle=\frac{1}{n_{i}}\boldsymbol{\phi}_{\alpha_{i}}^{\prime}(\boldsymbol{\Delta}_{i}\boldsymbol{\phi}_{\mu}\boldsymbol{\mu}_{k_{i}}+\boldsymbol{\Delta}_{i}\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k_{i}}+\boldsymbol{\phi}_{\alpha_{i}}\boldsymbol{\alpha}_{i}+\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i}-\boldsymbol{Y}_{i})=\boldsymbol{0},$
$\displaystyle\frac{\partial\;\textsf{DPSS}}{\partial\;\boldsymbol{\gamma}_{i}}$
$\displaystyle=\frac{1}{n_{i}}\boldsymbol{R}_{U_{i1}}(\boldsymbol{\Delta}_{i}\boldsymbol{\phi}_{\mu}\boldsymbol{\mu}_{k_{i}}+\boldsymbol{\Delta}_{i}\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k_{i}}+\boldsymbol{\phi}_{\alpha_{i}}\boldsymbol{\alpha}_{i}+\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i}-\boldsymbol{Y}_{i})+\lambda_{U_{i}}\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i}=\boldsymbol{0},$

which lead to

$\displaystyle\boldsymbol{\phi}_{\mu}^{\prime}\boldsymbol{\Delta}\boldsymbol{\phi}_{\mu}\boldsymbol{\mu}_{k}+\boldsymbol{\phi}_{\mu}^{\prime}\boldsymbol{\Delta}\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k}$
$\displaystyle=\boldsymbol{\phi}_{\mu}^{\prime}\boldsymbol{\tilde{Y}}_{k},$
(6)
$\displaystyle\boldsymbol{R}_{M_{1}}\boldsymbol{\Delta}\boldsymbol{\phi}_{\mu}\boldsymbol{\mu}_{k}+(\boldsymbol{R}_{M_{1}}\boldsymbol{\Delta}+\lambda_{M_{k}}\boldsymbol{I})\boldsymbol{R}_{M_{1}}\boldsymbol{\nu}_{k}$
$\displaystyle=\boldsymbol{R}_{M_{1}}\boldsymbol{\tilde{Y}}_{k},$ (7)
$\displaystyle\boldsymbol{\phi}_{\alpha_{i}}^{\prime}\boldsymbol{\phi}_{\alpha_{i}}\boldsymbol{\alpha}_{i}+\boldsymbol{\phi}_{\alpha_{i}}^{\prime}\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i}$
$\displaystyle=\boldsymbol{\phi}_{\alpha_{i}}^{\prime}\boldsymbol{\tilde{Y}}_{i},$
(8)
$\displaystyle\boldsymbol{R}_{U_{i1}}\boldsymbol{\phi}_{\alpha_{i}}\boldsymbol{\alpha}_{i}+(\boldsymbol{R}_{U_{i1}}+n_{i}\lambda_{U_{i}}\boldsymbol{I})\boldsymbol{R}_{U_{i1}}\boldsymbol{\gamma}_{i}$
$\displaystyle=\boldsymbol{R}_{U_{i1}}\boldsymbol{\tilde{Y}}_{i},$ (9)

with
$\boldsymbol{\tilde{Y}}_{k}=\displaystyle\sum_{i:k_{i}=k}\frac{1}{n_{i}}\boldsymbol{\Delta}_{i}^{\prime}\left(\boldsymbol{Y}_{i}-\boldsymbol{\phi}_{\alpha_{i}}\boldsymbol{{\alpha}}_{i}-\boldsymbol{R}_{U_{i1}}\boldsymbol{{\gamma}}_{i}\right)$,
$\boldsymbol{\tilde{Y}}_{i}=\boldsymbol{Y}_{i}-\boldsymbol{\Delta}_{i}\boldsymbol{\phi}_{\mu}\boldsymbol{{\mu}}_{k_{i}}-\boldsymbol{\Delta}_{i}\boldsymbol{R}_{M_{1}}\boldsymbol{{\nu}}_{k_{i}}$
and
$\boldsymbol{\Delta}=\displaystyle\sum_{i:k_{i}=k}\frac{1}{n_{i}}\boldsymbol{\Delta}_{i}^{\prime}\boldsymbol{\Delta}_{i}$.
The solutions of $\boldsymbol{\alpha}_{i}$ and $\boldsymbol{\gamma}_{i}$ in
the step (a) can be obtained from equations (8) and (9), while the solutions
of $\boldsymbol{\mu}_{k}$ and $\boldsymbol{\nu}_{k}$ in the step (b) from
equations (6) and (7).

###  A.3 Proof of Proposition 1

Based on the SDE $D^{p}X(t)=\dot{W}(t)$, we have

$D^{1}\boldsymbol{X}(t)=\boldsymbol{C}\boldsymbol{X}(t)+\boldsymbol{D}\dot{W}(t),$

where $\boldsymbol{X}=\\{X(t),D^{1}X(t),\cdots,D^{p-1}X(t)\\}^{\prime}$,
$\boldsymbol{C}=(c_{ii^{\prime}})_{p\times p}$, $c_{ii^{\prime}}=1$ when
$i^{\prime}=i+1$ and $c_{ii^{\prime}}=0$ otherwise, and
$\boldsymbol{D}=(0,0,\cdots,1)^{\prime}$.

It follows that

$\displaystyle\boldsymbol{X}_{j+1}$
$\displaystyle=\exp(\boldsymbol{C}\delta_{j})\boldsymbol{X}_{j}+\int_{0}^{\delta_{j}}\exp\\{\boldsymbol{C}(\delta_{j}-u)\\}\boldsymbol{D}\boldsymbol{\dot{W}}(t_{j}+u)du$
$\displaystyle=\boldsymbol{G}_{j}\boldsymbol{X}_{j}+\boldsymbol{\omega}_{j},$

where
$\boldsymbol{G}_{j}=\exp(\boldsymbol{C}\delta_{j})=\sum_{k=0}^{p}\frac{\delta_{j}^{k}}{k!}\boldsymbol{C}^{k}$
and
$\boldsymbol{\omega}_{j}\sim\textsf{N}_{p}\left(\boldsymbol{0},\boldsymbol{W}_{j}\right)$
with

$\displaystyle\boldsymbol{W}_{j}$
$\displaystyle=\int_{0}^{\delta_{j}}\exp\\{\boldsymbol{C}(\delta_{j}-u)\\}\boldsymbol{D}\boldsymbol{D}^{\prime}\exp\\{\boldsymbol{C}^{\prime}(\delta_{j}-u)\\}du$

as required.

## References

  * Barndorff-Nielsen and Shephard (2012) Barndorff-Nielsen, O. and Shephard, N. (2012), Financial volatility in continuous time, Cambridge: Cambridge University Press.
  * Durante et al. (2012) Durante, D., Scarpa, B., and Dunson, D. (2012), “Locally adaptive Bayesian covariance regression,” Arxiv preprint, arXiv:1210.2022.
  * Durbin and Koopman (2002) Durbin, J. and Koopman, S. (2002), “A simple and efficient simulation smoother for state space time series analysis,” Biometrika, 89, 603–616.
  * Durbin and Koopman (2001) Durbin, J. and Koopman, S.J. (2001), Time Series Analysis by State Space Methods, Oxford: Oxford University Press.
  * Heston (1993) Heston, S. (1993), “A closed-form solution for options with stochastic volatility with applications to bond and currency options,” Review of Financial Studies, 6, 327–343.
  * Ishihara and Omori (2012) Ishihara, T. and Omori, Y. (2012), “Efficient Bayesian estimation of a multivariate stochastic volatility model with cross leverage and heavy-tailed errors,” Computational Statistics & Data Analysis, 56, 3674–3689.
  * Jacquier et al. (2002) Jacquier, E., Polson, N., and Rossi, P. (2002), “Bayesian analysis of stochastic volatility models,” Journal of Business and Economic Statistics, 20, 69–87.
  * Jones (1993) Jones, R.H. (1993), Longitudinal Data with Serial Correlation: a State-space Approach, New York: Chapman & Hall/CRC.
  * Loddo et al. (2011) Loddo, A., Ni, S., and Sun, D. (2011), “Selection of Multivariate Stochastic Volatility Models via Bayesian Stochastic Search,” Journal of Business & Economic Statistics, 29, 342–355.
  * Miranda et al. (2009) Miranda, M.L., Maxson, P., and Edwards, S. (2009), “Environmental contributions to disparities in pregnancy outcomes,” Epidemiologic Reviews, 31, 67.
  * Müller and Yao (2010) Müller, H.G. and Yao, F. (2010), “Empirical dynamics for longitudinal data,” The Annals of Statistics, 38, 3458–3486.
  * Park and Choi (2010) Park, S. and Choi, S. (2010), “Hierarchical Gaussian process regression,” in Asian Conference on Machine Learning, pp. 95–110.
  * Raimann et al. (2012) Raimann, J., Usvyat, L., Thijssen, S., Kotanko, P., Rogus, J., Lacson, E., and Levin, N. (2012), “Blood pressure stability in hemodialysis patients confers a survival advantage: results from a large retrospective cohort study,” Kidney International, 81, 548–558.
  * Shephard (2005) Shephard, N. (2005), Stochastic volatility: selected readings, Oxford: Oxford University Press.
  * Van Es and Spreij (2011) Van Es, B. and Spreij, P. (2011), “Estimation of a multivariate stochastic volatility density by kernel deconvolution,” Journal of Multivariate Analysis, 102, 683–697.
  * Wahba (1990) Wahba, G. (1990), Spline models for observational data, vol. 59, Philadelphia: Society for Industrial Mathematics.
  * Wang et al. (2008) Wang, S., Jank, W., Shmueli, G., and Smith, P. (2008), “Modeling price dynamics in eBay auctions using differential equations,” Journal of the American Statistical Association, 103, 1100–1118.
  * Yao et al. (2005) Yao, F., Müller, H., and Wang, J. (2005), “Functional linear regression analysis for longitudinal data,” The Annals of Statistics, 33, 2873–2903.
  * Zhu et al. (2011) Zhu, B., Taylor, J., and Song, P. (2011), “Semiparametric Stochastic Modeling of the Rate Function in Longitudinal Studies,” Journal of the American Statistical Association, 106, 1485–1495.