{'fraction_non_alphanumeric': 0.04289322074027307, 'fraction_numerical': 0.014718648883276774, 'mean_word_length': 4.481938135409907, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stochastic Volatility Regression for Functional Data Dynamics

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.


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 1(a) 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 1(b). 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.
[ Figure 1 about here.] Additionally, empirical volatilities in Figure 1(b) 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.

The Model Specification
Suppose that Y i (t), i = 1, 2, . . . , m, is the observation of the ith subject at time t ∈ T i = {t i,1 , t i,2 , · · · , t i,n i < t U } with T i the set of observation times before time t U for the ith subject.
We specify an observation equation for Y i (t) as where Y i (t) is contaminated by the measurement error ε i (t) following a one-dimensional normal distribution with mean 0 and variance σ 2 ε . Assuming the ith subject belongs to the k i th group (e.g. by race or treatment) with k i ∈ {1, 2, · · · , g}, we include a k i th group mean function M k i (t) = E{Y i (t) | M k i (t)} in the observation equation. In addition, the trajectory of the ith 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 E{U i (t)} = 0.
The volatility of the ith 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 σ 2 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 σ 2 M k i and σ 2 U i : where p, q ∈ N ≥ 1 and and U i (t) as well as their derivatives till orders q − 1 and p − 1 respectively. The volatility σ 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(σ 2 U i ) ∼ N 1 (x ′ i β, σ 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, . . . , g, and U i (t), i = 1, 2, . . . , n, are the summations of mutually independent Gaussian processes written as Hence, we can represent the prior of M k i (t) + U i (t) as a hierarchical Gaussian process, where GP(M(t), K(s, t)) denotes a Gaussian process with mean function M(t) and covariance function 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 is subject specific and depends on covariates through σ 2 U i . To carry out Bayesian inference, we further specify the following prior distributions for the hyperparameters. In particular, ) denotes the inverse gamma distribution with shape parameter a and scale parameter b. The β and σ 2 follow the independent Jeffreys' prior, f (β, σ 2 ) ∝ 1/σ 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 σ 2 M 0 and σ 2 U 0 go to infinity, the posterior means of M k (t) and U i (t) are equivalent to the double-penalized smoothing splineM k (t) +Û i (t), which is the minimizer of the double-penalized sum-of-squares, U i (t) respectively, where the smoothness and the fidelity to data are balanced by the smoothing . Expression forM k (t) andÛ i (t) can be obtained explicitly, as shown in the following Theorem.
Theorem 1 The smoothing splinesM k (t) andÛ i (t) with t ∈ T minimize the double-penalized sum-of-squares (4) and have the formŝ . . , n}, the index set of unique observation times among all m subjects.
The proofs of Theorem 1 and the following Corollary are included in Appendix A.
Corollary 1 The µ k , ν k , α i and γ i can be obtained through a backfitting algorithm or the Gauss-Seidel method, iterating the following two steps until convergence:

Posterior Computation
Although we can obtainM k (t) andÛ i (t) by the backfitting algorithm outlined in Corollary 1 and estimate λ M k and λ U i through generalized cross validation (Chap. 4, Wahba, 1990), it is unclear how to allow λ U i to depend on covariates. In addition, when n is large, it is computational infeasible to invert the n × n matrix 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) =Ẇ (t). Consequently, the X j = {X(t j ), D 1 X(t j ), · · · , D r−1 X(t j )} ′ , j = 1, 2, · · · , n, follows a state equation The proof is in Appendix A. Finally, we outline the proposed MCMC as follows.
(1) Given M k i (t ij ), σ 2 ε and σ 2 U i , sample U i (t ij ), i = 1, 2, · · · , m, j = 0, 1, · · · , n i . Let Y U ij = Y i (t ij ) − M k i (t ij ) and the SVR model for the ith 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.
. Similar to the G j , ω j and W j in Proposition 1, the G U ij , ω U ij and W U ij follow the same specifications with r = q.
(2) Given U i (t j ), σ 2 ε and σ 2 M k , sample M k (t j ), k = 1, 2, · · · , g, j = 0, 1, · · · , n. Similarly, we rewrite the SVR model for the kth group as the following state space model and then sample M k i (t ij ) and its derivatives by the simulation smoother.
When ith subject has an observation at time t j and k i = k, The G M kj , ω M kj and W M kj are given by Proposition 1 with r = p.
(3a) Given M k i (t ij ) and U i (t ij ), i = 1, 2, · · · , m, j = 1, 2, · · · , n i , sample σ 2 ε from the posterior distribution invGamma a + 1 (3d) Given U ij , β and σ 2 , sample σ 2 U i using a Metropolis-Hasting algorithm. We choose σ 2 U i ∼ invGamma(a, b) as the proposal prior distribution and a proposal σ 2 * U i can be easily drawn from invGamma a + n i q the corresponding proposal posterior distribution. The σ 2 * U i will be accepted with the following probability and discarded otherwise with σ 2 U i unchanged, where f LN , f Nq and f iG denote the log-normal, q-dimensional normal and inverse gamma probability density functions respectively with a U i = a+ n i q We draw τ from Chi-squared distribution with m − k degrees of freedom and set σ 2 = (m−k)σ 2 τ and then sample β from N m β , σ 2 (X ′ X) −1 .

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 σ 2 U i and the coefficients β. 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 β = (0, 0.6, 2) ′ and ∼ N 1 (0, 0.25) respectively. Given β and x i , volatilities σ 2 U i 's can be drawn from log(σ 2 U i ) ∼ N 1 (x ′ i β, 1). Along with σ 2 M 1 = σ 2 M 2 = 10, σ 2 ε = 1, p = 2 and q = 1, M 1 (t), M 2 (t), U i (t) and ε i (t) are sampled at t ∈ {0.2, 0.4, · · · , 4} from equations (2) and (3) and the distribution of measurement error ε i (t).
Twenty percent of samples are removed completely at random, resulting in on average 16 unequally spaced observations per subject. Finally, the ith 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), σ 2 U i and β. 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 1 in the first stage withÛ 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 β.
For each simulated dataset, we calculate average squared error (ASE) for the trajectory ASE(M + , and squared errors (SE) for coefficient estimates SE(β l ) = (β l −β l ) 2 , l = 0, 1, 2. Table 1 reports means of ASE(M + U), ASE{log(σ 2 U )} and SE(β l ) across 100 replicate datasets. MASEs and MSEs by NCS and FPCA approaches are significantly inflated, for example, being doubled and tripled in 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 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 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 ∈ {0.2, 0.4, · · · , 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α 1i cos(πt/10) + 0.2α 2i sin(πt/10) + ε i (t) for subjects in the first group and from Y i (t) = 10{t + cos(t)} + 0.5α 1i cos(πt/10) + 0.3α 2i sin(πt/10) + ε i (t) for the ones in the second group, with α 1i ∼ N 1 (0, 1). The SVR, NCS and FPCA approaches are applied and MASE(M + U) is presented in Table 1, in which SVR and FPCA approaches show close 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.

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 ≤ 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 (a) and (b) 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 (c) 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 (d). 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.

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.
Hence, the double-penalized sum-of-squares (4) can be written as

A.2 Proof of Corollary 1
Taking partial derivatives of double penalized sum-of-squares in Corollary 1 with respective to µ k , ν k , α i and γ i respectively and setting them to zeros, we have The solutions of α i and γ i in the step (a) can be obtained from equations (8) and (9), while the solutions of µ k and ν k in the step (b) from equations (6) and (7).
It follows that as required.  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 ith woman at time t ij , and U ij = Y ij −Ȳ j is the deviation from the group mean blood pressureȲ j . The empirical volatility measures the fluctuation of trajectories empirically and is defined as 1 with n i the number of observations for the ith woman. The plots of observation (•) and trajectory at time t ij (×), as well as estimates of trajectory M k i (t) + U i (t) by SVR (-), NCS (− − −) and FPCA (· · ·) approaches, for six subjects in one simulated dataset with the largest individual average squared errors