Modeling Heart Rate Variability (HRV) data has become important for clinical applications and as a research tool. These data exhibit long memory and time-varying conditional variance (volatility). In HRV, volatility is traditionally estimated by recursive least squares combined with short memory AutoRegressive (AR) models. This work considers a parametric approach based on long memory Fractionally Integrated AutoRegressive Moving Average (ARFIMA) models with heteroscedastic errors. To model the heteroscedasticity nonlinear Generalized Autoregressive Conditionally Heteroscedastic (GARCH) and Exponential Generalized Autoregressive Conditionally Heteroscedastic (EGARCH) models are considered. The latter are necessary to model empirical characteristics of conditional volatility such as clustering and asymmetry in the response, usually called leverage in time series literature. The ARFIMA-EGARCH models are used to capture and remove long memory and characterize conditional volatility in 24 hour HRV recordings from the Noltisalis database. © 2016 IEEE

Leite, A

Maria Eduarda Silva

Rocha, AP

Repositório Aberto da Universidade do Porto

Modeling volatility in Heat Rate VariabilityArgentina Leite1 and Maria Eduarda Silva2 and Ana Paula Rocha3,Member, IEEEAbstract—Modeling Heart Rate Variability (HRV) data hasbecome important for clinical applications and as a researchtool. These data exhibit long memory and time-varying condi-tional variance (volatility). In HRV, volatility is traditionallyestimated by recursive least squares combined with shortmemory AutoRegressive (AR) models. This work considers aparametric approach based on long memory Fractionally Inte-grated AutoRegressive Moving Average (ARFIMA) models withheteroscedastic errors. To model the heteroscedasticity nonlin-ear Generalized Autoregressive Conditionally Heteroscedastic(GARCH) and Exponential Generalized Autoregressive Con-ditionally Heteroscedastic (EGARCH) models are considered.The latter are necessary to model empirical characteristics ofconditional volatility such as clustering and asymmetry in theresponse, usually called leverage in time series literature. TheARFIMA-EGARCH models are used to capture and removelong memory and characterize conditional volatility in 24 hourHRV recordings from the Noltisalis database.I. INTRODUCTIONHeart Rate Variability (HRV) reflects the interaction be-tween perturbations to the cardiovascular variables and thecorresponding response of the cardiovascular regulatory sys-tems [1]. Thus modeling such variability can provide a quan-titative and non-invasive method to assess the integrity ofthe cardiovascular system. HRV data display non stationarycharacteristics and exhibit long memory and time-varyingconditional variance (usually designated by volatility) whichmay contain indicators of current disease or warnings aboutimpending diseases.Traditionally, HRV data can be characterized by linearAutoRegressive (AR) models, which describe only shortmemory in the mean. These models combined with recur-sive least squares have been used to estimate volatility inHRV data [2]. However, it is acknowledged that complexinteractions of electrophysiological, humoral variables andautonomic and central nervous regulations induce nonlineareffects in HRV [3], [4]. Thus the analysis of HRV basedon nonlinear models might elicit valuable information for*This work is financed by the ERDF - European Regional Develop-ment Fund through the Operational Programme for Competitiveness andInternationalisation - COMPETE 2020 Programme within project POCI-01-0145-FEDER-006961, and by National Funds through the Portuguesefunding agency, FCT - Fundac¸a˜o para a Cieˆncia e a Tecnologia as part ofprojects UID/EEA/50014/2013, CIDMA UID/MAT/04106/2013 and CMUPUID/MAT/00144/2013, funded by FCT (Portugal) with national (MEC)and European structural funds through the programs FEDER, under thepartnership agreement PT2020.1Argentina Leite is with Escola de Cieˆncias e Tecnologia, Universidadede Tra´s-os-Montes e Alto Douro & C-BER & INESC TEC, Portugaltinucha@utad.pt2Maria Eduarda Silva is with Faculdade de Economia, Universidade doPorto & CIDMA, Portugal mesilva@fep.up.pt3Ana Paula Rocha is with Faculdade de Cieˆncias, Universidade do Porto& CMUP, Portugal aprocha@fc.up.ptits physiological interpretation as well as risk assessment.A nonlinear approach to volatility description in HRV datawas proposed by Leite et al [5], [6] using Fractionally Inte-grated AutoRegressive Moving Average (ARFIMA) modelswith Generalized AutoRegressive Conditional Heteroscedas-tic (GARCH) errors, an extension of the usual AR anal-ysis. ARFIMA-GARCH models are used to capture andremove long memory and estimate volatility in 24 hour HRVrecordings [6], [7]. However, GARCH models assume thatvolatility depends only on the magnitude of the shocks andnot on their sign, meaning that positive and negative shockshave a symmetric effect on volatility [8].This work considers Exponential GARCH (EGARCH)[8] models which are an extension of GARCH models.EGARCH models assume that the effect of positive andnegative shocks on volatility is asymmetric, an effect usuallydesignated by leverage. A model with leverage effect is moresuited to describe the complex characteristics of HRV data.The ARFIMA-EGARCH models are applied to 24 hour HRVrecordings from 30 subjects of the Noltisalis database [9]: 10healthy subjects (N), 10 patients suffering from congestiveheart failure (C) and 10 heart transplanted patients (T). Thedatabase was collected by the cooperative effort of universitydepartments and rehabilitation clinics in Italy, [9].II. MODELS WITH CONDITIONAL VOLATILITYIn this work we consider models satisfyingφ(B)(1−B)dxt = θ(B)t (1)t = σtzt (2)σ2t = Var(t|Ht−1) (3)where B is the backward-shift operator. Equation (1) de-scribes the conditional mean of the process with seriallyuncorrelated residuals t and is said an ARFIMA(p, d, q)with p, q ∈ N0 and d ∈ R, [10]; d is the long-memoryparameter, determines the long-term behaviour in mean and(1−B)d =∞∑k=0(dk)(−1)kBk is the fractional differenceoperator [10]; p, q ∈ N and the polynomials φ(B) = 1 −φ1B − ...− φpBp and θ(B) = 1 + θ1B + ...+ θqBq allowfor the modeling of the short-range properties in the mean;for −0.5 < d < 0.5 the ARFIMA process is covariancestationary. In the range −0.5 < d < 0.5, the long memoryparameter is related to the Hurst coefficient, H , to the fractaldimension, D, and to the slope of the (generalized) spectraldensity in the low frequency range, α, by d = H−0.5, H =2−D and α = 2d, respectively. Moreover, for 0.5 ≤ d < 1the process ARFIMA is non-stationary and mean reverting.Equations (2) and (3) describe the conditional variance ofthe process which varies over time as in time-varying ARmodels. In (2) t are called shocks and zt, independent andidentically distributed random variables with zero mean andunit variance, are the standardized shocks. The conditionalvariance σ2t in (2) depends on Ht−1 which is the historyof process and includes past variances σ21 , σ22 , . . . , σ2t−1, andpast shocks 1, 2, . . . , t−1.There are several models to govern the evolution of σ2t .The most common is the GARCH(P,Q), P,Q ∈ N0 model[8] under whichσ2t = u0 +P∑i=1vi σ2t−i +Q∑i=1ui 2t−i (4)This model requires non-negativity constraints on the pa-rameters to ensure positive conditional variances, u0 >0, v1, ..., vP , u1, ..., uQ ≥ 0,P∑i=1vi +Q∑j=1uj < 1. Theparameters, ui and vi characterize the volatility clusteringphenomena observed in many data sets. It is noteworthy thatGARCH models assume that conditional volatility dependsonly on magnitude of the shocks and not on their sign.Several extensions of GARCH models have been proposedso that empirical characteristics of the data such as asymmet-ric response to positive and negative shocks, are adequatelymodeled. One such model is the EGARCH(P,Q) [8] definedas follows:logσ2t = u∗0+P∑i=1vi logσ2t−i+Q∑i=1ui |zt−i|+Q∑i=1ξizt−i (5)where u∗0 = u0 −Q∑i=1ui√2pi and zt = t/σt.This process does not require constraints on the parametersfor ensuring the positivity of the variance. The parameters uiand vi characterise the volatility clustering phenomena andthe parameters ξi describe the leverage effects. The impactof positive shocks, t−i > 0 is (ui + ξi)t−iσt−i, while fornegative shocks it is (ui−ξi) t−iσt−i . If ξi = 0, logσ2t respondssymmetrically to t−i.In this work, we consider ARFIMA(p, d, 0)-GARCH(1, 1)and ARFIMA(p, d, 0)-EGARCH(1, 1) models, since they area natural extension of the classic AR(p) models usual in theanalysis of HRV. Additionally to the parameters d, whichcharacterises the long memory in the mean, u1 and v1,which characterize the volatility clustering, we give specialattention to ξ which describes the leverage effect in theconditional variance.Given a time series x1, . . . , xn, to estimate the parametersof the above mentioned models proceed as follows [8],[10]: (i) estimate d using the semi-parametric local Whittleestimator; (ii) define the filtered data yt = (1−B)dxt; (iii)estimate the AR(p)-GARCH(1, 1) or AR(p)-EGARCH(1, 1)parameters in the filtered data yt by maximum likelihood(Econometrics Toolbox of MATLAB [11]), with the order pdetermined by the Akaike Information Criterion (AIC).Fig. 1. Short-term HRV data: (a) tachogram of a normal subject (segmentwith 1700 beats of RR series for subject-N6 from the Noltisalis database),(b) ACF of the data, (c) residuals of the fitted ARFIMA(9,0.34,0) model,(d) ACF of the residuals, (e) ACF of the squared residuals and (f) CCF ofthe residuals and the squared residuals. The horizontal lines (- -) show the95% confidence limits, (g) estimated volatility, σˆt.III. MODELING HRVTo motivate the use of ARFIMA(p,d,0)-EGARCH(1,1)models in HRV data, consider the tacogram of the healthysubject N6 (segment of RR series with 1700 beats) rep-resented in Fig. 1(a). The ACF (autocorrelation function)of the data, represented in (b), shows a typical very slowdecay indicating the presence of long memory which canbe adequately modelled with an ARFIMA, equation (1).The residuals (ˆt) from the fitted ARFIMA(p = 9, d =0.34, 0) (chosen by AIC) are represented in Fig. 1(c). Thecorresponding ACF in (d), exhibits small correlation indi-cating that the ARFIMA model is adequate to explain thedynamics of the conditional mean of the data. However, thesquared residuals exhibit significant autocorrelation in Fig.1(e), indicating time-varying conditional variance which canbe modeled by a GARCH model, (4), hereafter denoted asModel I: ARFIMA(9,0.34,0)-GARCH(1,1). Additionally, thecross-correlation between the residuals ˆt and their squaresˆ2t , Fig. 1(f), indicates that the conditional variance dependsalso on the lagged shocks, an effect that can be modeledby an EGARCH model, (5), hereafter denoted as ModelII: ARFIMA(9,0.34,0)-EGARCH(1,1). The results for twomodels are summarized in Table I.McLeod-Li testing [10] of the residuals of Model II gives atest statistic Q(20) = 41.40, indicating no significant condi-tional heteroscedasticity in the residual series. Moreover, thecross-correlation between the residuals and their squares inlag 1, CCF (1), is 0.04 indicating no significant correlationand the AIC criterion favours Model II.These results indicate that the ARFIMA(9,d,0)-EGARCH(1,1) (Model II) model leads to furthercharacterisation of the HRV data. In fact, the estimateξˆ1 = 0.24 indicates asymmetric response to positiveand negative shocks. The conditional standard deviationestimate σˆt represented in Fig. 1(g), captures very well theheteroscedasticity in the original data, plotted in (a). Similarresults were obtained in other short HRV recordings.TABLE IESTIMATES (STANDARD DEVIATION) FOR MODEL I AND MODEL IIFITTED TO DATA IN FIG. 1. Q(20) STANDS FOR MCLEOD-LI’ S TESTSTATISTIC FOR THE RESIDUALS.Parameter Model I Model IIdˆ 0.34 0.34Q(20) of ˆt 494.32 494.32CCF(1) of ˆt,ˆ2t 0.35 0.35uˆ0 0.13e−03(0.16e−04) −2.87 (0.20)uˆ1 0.33 (0.02) 0.47 (0.05)vˆ1 0.40 (0.05) 0.63 (0.03)ξˆ1 —– 0.24 (0.02)Q(20) of zˆt 28.84 41.40CCF(1) of zˆt,zˆ2t 0.12 0.04AIC of zˆt -8578.00 -8623.00In the case of long recordings, such as ambulatory 24hour HRV data (approximately 100000 beats), exhibitingseveral non stationary characteristics, ARFIMA-EGARCHmodeling is combined with adaptive segmentation [6]: longrecords are decomposed into short records of variable lengthand the break points, which mark the end of consecutiveshort records, are identified by AIC criterion. The shortrecords thus obtained have a minimum length 512 and aresubsequently modeled using ARFIMA-EGARCH models.IV. RESULTS AND DISCUSSIONIn this section, the above methodology is applied to long-term HRV series of subjects from the Noltisalis database [9].The results are first illustrated for the healthy subject N6 andthe patient C10, in Figs. 2 and 3, respectively.RR(sec)0.511.5(a)d00.51(b)Q(20)0200400600(c)CCF(1)00.20.4(d)u1, v1(--)00.51(e)Recording time (hh:mm)10:00 15:00 20:00 01:00 06:00 11:00ξ100.20.40.6(f)Fig. 2. (a) Tachogram of healthy subject N6, 24 hour recordings Noltisalisdatabase, Evolution over 24 hours of dˆ in (b), uˆ1(-) and vˆ1(- -) in (e) andξˆ1 in (f) estimated using ARFIMA-EGARCH models and segmentation;statistic Q(20) of the McLeod-Li test for conditional heteroscedasticity in(c) and cross-correlation between the residuals t and their squares of lag1 in (d).The long memory estimates dˆ in Fig. 2(b) and Fig. 3(b),change over time showing circadian variation, with lowestvalues during the night periods, 0 < dˆ < 0.5 in contrastwith 0.5 < dˆ < 1 for the day period. These results are inconcordance with Leite et al.[6].The residuals from the ARFIMA modeling are tested forconditional heteroscedasticity with McLeod-Li test. The cor-responding test statistics Q(20) are represented in plots (c)in Figs. 2 and 3. It is found that the percentage of segmentswith heteroscedasticity are 98.8% and 94.5 % for subjectsN6 and C10, respectively. Moreover, the volatility parametersestimates uˆ1 and vˆ1, Fig. 2(e) and Fig. 3(e), change over timewith some circadian variation for the healthy subject. Forthe sick subject, the estimate vˆ1 increases and uˆ1 decreases.Note that the estimated values for parameter v1 are over 0.5indicating some persistence in conditional variance. Theseresults are in agreement with Leite et al.[6].The residuals from the ARFIMA models are also checkedfor the asymmetric effect (leverage). Plot (d) in Figs. 2 and3 represent the values of CCF(1) between the residuals andtheir squares. The correlations indicate that the conditionalvariance depends also on the lagged shocks, leverage effect.In fact, the percentage of segments with leverage effect are93.8% and 43.8 % for subjects N6 and C10, respectively.RR(sec)0.511.5(a)d00.51 (b)Q(20)0200400600 (c)CCF(1)00.20.4(d)u1, v1(--)00.51(e)Recording time (hh:mm)10:00 15:00 20:00 01:00 06:00 11:00ξ100.20.4(f)Fig. 3. (a) Tachogram of patient C10 affected by congestive heart failure, 24hour recordings Noltisalis database. Evolution over 24 hours of dˆ in (b), uˆ1(-) and vˆ1(- -) in (e) and ξˆ1 in (f) estimated using ARFIMA-EGARCH modelsand segmentation; statistic Q(20) of the McLeod-Li test for conditionalheteroscedasticity in (c) and cross-correlation between the residuals t andtheir squares of lag 1 in (d).Moreover, the leverage parameter estimates ξˆ1, Fig. 2(f) andFig. 3(f), change over time and present higher values for thehealthy subject. The same analysis is performed for patientT3. The results for subjects N6, C10 and T3 are summarizedin Fig. 4(a).Finally, the mean estimates dˆ, uˆ1, vˆ1 and ξˆ1 for each ofthe 30 subjects in the database are summarized in Fig. 4(b).These overall results indicate that the long memory param-eter d is lower for healthy subjects N while the volatilityand the leverage parameters u1, ξ1 respectively, are lowerfor sick subjects C and T. The parameter v1 presents highvariability for the sick subjects T. The results are promisingin differentiating health and disfunction situations deservingfurther study.V. CONCLUSIONSThis paper contributes for further characterization of thecomplex dynamics of HRV with a leverage parameter. Inparticular, the model under discussion indicates that valuesof HRV under the mean lead to less variability in HRV thanvalues over the mean- asymmetric effect captured by theleverage parameter. Furthermore this study indicates that thiseffect is stronger in healthy subjects. However, further studiesare necessary to assess the importance of this parameter inN6 C10 T3 N6 C10 T3 N6 C10 T3 N6 C10 T3Estimated-0.200.20.40.60.81(a)N C T N C T N C T N C TEstimated-0.200.20.40.60.81(b)d u1v1ξ1ξ1v1u1dFig. 4. Boxplots for dˆ, uˆ1, vˆ1 and ξˆ1 using ARFIMA-EGARCH modelsand segmentation: (a) an healthy subject-N6, a patient suffering fromcongestive heart failure-C10 and a heart transplanted patient-T3, (b) subjectsmean values for each group N, C, T.risk stratification and autonomic nervous system dysfunctioncharacterization.REFERENCES[1] Task Force of the European Society of Cardiology and North AmericanSociety of Pacing Electrophysiology, Heart rate variability: standardsof measurement, physiological interpretation and clinical use, Circu-lation, vol. 93, pp. 1043-1065, 1996.[2] L. T. Mainardi, On the quantification of heart rate variability spectralparameters using time-frequency and time-varying methods, Phil.Trans. R. Soc. A, vol. 367, pp. 2555-275, 2009.[3] R. T. Baillie, A. A. Cecen, and C. Erkal, Normal heartbeat series arenonchaotic, nonlinear, and multifractal: New evidence from semipara-metric and parametric tests, Chaos, vol. 19, pp. 028503, 2009.[4] A. Voss, S. Schulz, R. Schroeder, M. Baumert, and P. Caminal,Methods derived from nonlinear dynamics for analysing heart ratevariability. Phil. Trans. R. Soc. A, vol. 367, pp. 277-296, 2009.[5] A. Leite, A. P. Rocha, and M. E. Silva, Long memory and volatility inHRV: an ARFIMA-GARCH approach, in Computing in Cardiology,vol. 36, pp. 165-168, 2009.[6] A. Leite, A. P. Rocha, and M. E. Silva, Beyond long memory inheart rate variability: an approach based on fractionally integratedautoregressive moving average time series models with conditionalheteroscedasticity, Chaos, vol. 23, pp. 023103, 2013.[7] M. Y. Hu, and C. Tsoukalas, Conditional volatility properties of sleepdisordered breathing, Comput. Biol. Med., vol. 36, pp. 303-312, 2006.[8] R. S. Tsay, Analysis of Financial Time Series (2nd ed.). Wiley-Interscience, 2005.[9] M. G. Signorini, R. Sassi, S. Cerutti, Working on the NOLTISALISDatabase: Measurement of nonlinear properties in heart rate variabilitysignals. In: Proceedings of IEEE-EMBS International Conference,Istanbul, Turkey (IEEE, Piscataway), pp. 547-550, 2001.[10] P. Boukhan, G. Oppenheim, and M. S. Taqqu, Theory and applicationsof long-range dependence, Boston: Birkha¨user, 2003.[11] MATLAB, version 8.5 (R2015a), The MathWorks Inc., Natick, Mas-sachusetts, 2015.

Modeling volatility in heat rate variability

Argentina Leite

Ana Paula Rocha

Crossref

Open Repository of the University of Porto

Modeling volatility in Heat Rate Variability
Argentina Leite1 and Maria Eduarda Silva2 and Ana Paula Rocha3,Member, IEEE
Abstract—Modeling Heart Rate Variability (HRV) data has
become important for clinical applications and as a research
tool. These data exhibit long memory and time-varying condi-
tional variance (volatility). In HRV, volatility is traditionally
estimated by recursive least squares combined with short
memory AutoRegressive (AR) models. This work considers a
parametric approach based on long memory Fractionally Inte-
grated AutoRegressive Moving Average (ARFIMA) models with
heteroscedastic errors. To model the heteroscedasticity nonlin-
ear Generalized Autoregressive Conditionally Heteroscedastic
(GARCH) and Exponential Generalized Autoregressive Con-
ditionally Heteroscedastic (EGARCH) models are considered.
The latter are necessary to model empirical characteristics of
conditional volatility such as clustering and asymmetry in the
response, usually called leverage in time series literature. The
ARFIMA-EGARCH models are used to capture and remove
long memory and characterize conditional volatility in 24 hour
HRV recordings from the Noltisalis database.
I. INTRODUCTION
Heart Rate Variability (HRV) reflects the interaction be-
tween perturbations to the cardiovascular variables and the
corresponding response of the cardiovascular regulatory sys-
tems [1]. Thus modeling such variability can provide a quan-
titative and non-invasive method to assess the integrity of
the cardiovascular system. HRV data display non stationary
characteristics and exhibit long memory and time-varying
conditional variance (usually designated by volatility) which
may contain indicators of current disease or warnings about
impending diseases.
Traditionally, HRV data can be characterized by linear
AutoRegressive (AR) models, which describe only short
memory in the mean. These models combined with recur-
sive least squares have been used to estimate volatility in
HRV data [2]. However, it is acknowledged that complex
interactions of electrophysiological, humoral variables and
autonomic and central nervous regulations induce nonlinear
effects in HRV [3], [4]. Thus the analysis of HRV based
on nonlinear models might elicit valuable information for
*This work is financed by the ERDF - European Regional Develop-
ment Fund through the Operational Programme for Competitiveness and
Internationalisation - COMPETE 2020 Programme within project POCI-
01-0145-FEDER-006961, and by National Funds through the Portuguese
funding agency, FCT - Fundac¸a˜o para a Cieˆncia e a Tecnologia as part of
projects UID/EEA/50014/2013, CIDMA UID/MAT/04106/2013 and CMUP
UID/MAT/00144/2013, funded by FCT (Portugal) with national (MEC)
and European structural funds through the programs FEDER, under the
partnership agreement PT2020.
1Argentina Leite is with Escola de Cieˆncias e Tecnologia, Universidade
de Tra´s-os-Montes e Alto Douro & C-BER & INESC TEC, Portugal
tinucha@utad.pt
2Maria Eduarda Silva is with Faculdade de Economia, Universidade do
Porto & CIDMA, Portugal mesilva@fep.up.pt
3Ana Paula Rocha is with Faculdade de Cieˆncias, Universidade do Porto
& CMUP, Portugal aprocha@fc.up.pt
its physiological interpretation as well as risk assessment.
A nonlinear approach to volatility description in HRV data
was proposed by Leite et al [5], [6] using Fractionally Inte-
grated AutoRegressive Moving Average (ARFIMA) models
with Generalized AutoRegressive Conditional Heteroscedas-
tic (GARCH) errors, an extension of the usual AR anal-
ysis. ARFIMA-GARCH models are used to capture and
remove long memory and estimate volatility in 24 hour HRV
recordings [6], [7]. However, GARCH models assume that
volatility depends only on the magnitude of the shocks and
not on their sign, meaning that positive and negative shocks
have a symmetric effect on volatility [8].
This work considers Exponential GARCH (EGARCH)
[8] models which are an extension of GARCH models.
EGARCH models assume that the effect of positive and
negative shocks on volatility is asymmetric, an effect usually
designated by leverage. A model with leverage effect is more
suited to describe the complex characteristics of HRV data.
The ARFIMA-EGARCH models are applied to 24 hour HRV
recordings from 30 subjects of the Noltisalis database [9]: 10
healthy subjects (N), 10 patients suffering from congestive
heart failure (C) and 10 heart transplanted patients (T). The
database was collected by the cooperative effort of university
departments and rehabilitation clinics in Italy, [9].
II. MODELS WITH CONDITIONAL VOLATILITY
In this work we consider models satisfying
φ(B)(1−B)dxt = θ(B)t (1)
t = σtzt (2)
σ2t = Var(t|Ht−1) (3)
where B is the backward-shift operator. Equation (1) de-
scribes the conditional mean of the process with serially
uncorrelated residuals t and is said an ARFIMA(p, d, q)
with p, q ∈ N0 and d ∈ R, [10]; d is the long-memory
parameter, determines the long-term behaviour in mean and
(1−B)d =
∞∑
k=0
(
d
k
)
(−1)kBk is the fractional difference
operator [10]; p, q ∈ N and the polynomials φ(B) = 1 −
φ1B − ...− φpBp and θ(B) = 1 + θ1B + ...+ θqBq allow
for the modeling of the short-range properties in the mean;
for −0.5 < d < 0.5 the ARFIMA process is covariance
stationary. In the range −0.5 < d < 0.5, the long memory
parameter is related to the Hurst coefficient, H , to the fractal
dimension, D, and to the slope of the (generalized) spectral
density in the low frequency range, α, by d = H−0.5, H =
2−D and α = 2d, respectively. Moreover, for 0.5 ≤ d < 1
the process ARFIMA is non-stationary and mean reverting.
Equations (2) and (3) describe the conditional variance of
the process which varies over time as in time-varying AR
models. In (2) t are called shocks and zt, independent and
identically distributed random variables with zero mean and
unit variance, are the standardized shocks. The conditional
variance σ2t in (2) depends on Ht−1 which is the history
of process and includes past variances σ21 , σ
2
2 , . . . , σ
2
t−1, and
past shocks 1, 2, . . . , t−1.
There are several models to govern the evolution of σ2t .
The most common is the GARCH(P,Q), P,Q ∈ N0 model
[8] under which
σ2t = u0 +
P∑
i=1
vi σ
2
t−i +
Q∑
i=1
ui 
2
t−i (4)
This model requires non-negativity constraints on the pa-
rameters to ensure positive conditional variances, u0 >
0, v1, ..., vP , u1, ..., uQ ≥ 0,
P∑
i=1
vi +
Q∑
j=1
uj < 1. The
parameters, ui and vi characterize the volatility clustering
phenomena observed in many data sets. It is noteworthy that
GARCH models assume that conditional volatility depends
only on magnitude of the shocks and not on their sign.
Several extensions of GARCH models have been proposed
so that empirical characteristics of the data such as asymmet-
ric response to positive and negative shocks, are adequately
modeled. One such model is the EGARCH(P,Q) [8] defined
as follows:
logσ2t = u
∗
0+
P∑
i=1
vi logσ2t−i+
Q∑
i=1
ui |zt−i|+
Q∑
i=1
ξizt−i (5)
where u∗0 = u0 −
Q∑
i=1
ui
√
2
pi and zt = t/σt.
This process does not require constraints on the parameters
for ensuring the positivity of the variance. The parameters ui
and vi characterise the volatility clustering phenomena and
the parameters ξi describe the leverage effects. The impact
of positive shocks, t−i > 0 is (ui + ξi)
t−i
σt−i
, while for
negative shocks it is (ui−ξi) t−iσt−i . If ξi = 0, logσ2t responds
symmetrically to t−i.
In this work, we consider ARFIMA(p, d, 0)-GARCH(1, 1)
and ARFIMA(p, d, 0)-EGARCH(1, 1) models, since they are
a natural extension of the classic AR(p) models usual in the
analysis of HRV. Additionally to the parameters d, which
characterises the long memory in the mean, u1 and v1,
which characterize the volatility clustering, we give special
attention to ξ which describes the leverage effect in the
conditional variance.
Given a time series x1, . . . , xn, to estimate the parameters
of the above mentioned models proceed as follows [8],
[10]: (i) estimate d using the semi-parametric local Whittle
estimator; (ii) define the filtered data yt = (1−B)dxt; (iii)
estimate the AR(p)-GARCH(1, 1) or AR(p)-EGARCH(1, 1)
parameters in the filtered data yt by maximum likelihood
(Econometrics Toolbox of MATLAB [11]), with the order p
determined by the Akaike Information Criterion (AIC).
Fig. 1. Short-term HRV data: (a) tachogram of a normal subject (segment
with 1700 beats of RR series for subject-N6 from the Noltisalis database),
(b) ACF of the data, (c) residuals of the fitted ARFIMA(9,0.34,0) model,
(d) ACF of the residuals, (e) ACF of the squared residuals and (f) CCF of
the residuals and the squared residuals. The horizontal lines (- -) show the
95% confidence limits, (g) estimated volatility, σˆt.
III. MODELING HRV
To motivate the use of ARFIMA(p,d,0)-EGARCH(1,1)
models in HRV data, consider the tacogram of the healthy
subject N6 (segment of RR series with 1700 beats) rep-
resented in Fig. 1(a). The ACF (autocorrelation function)
of the data, represented in (b), shows a typical very slow
decay indicating the presence of long memory which can
be adequately modelled with an ARFIMA, equation (1).
The residuals (ˆt) from the fitted ARFIMA(p = 9, d =
0.34, 0) (chosen by AIC) are represented in Fig. 1(c). The
corresponding ACF in (d), exhibits small correlation indi-
cating that the ARFIMA model is adequate to explain the
dynamics of the conditional mean of the data. However, the
squared residuals exhibit significant autocorrelation in Fig.
1(e), indicating time-varying conditional variance which can
be modeled by a GARCH model, (4), hereafter denoted as
Model I: ARFIMA(9,0.34,0)-GARCH(1,1). Additionally, the
cross-correlation between the residuals ˆt and their squares
ˆ2t , Fig. 1(f), indicates that the conditional variance depends
also on the lagged shocks, an effect that can be modeled
by an EGARCH model, (5), hereafter denoted as Model
II: ARFIMA(9,0.34,0)-EGARCH(1,1). The results for two
models are summarized in Table I.
McLeod-Li testing [10] of the residuals of Model II gives a
test statistic Q(20) = 41.40, indicating no significant condi-
tional heteroscedasticity in the residual series. Moreover, the
cross-correlation between the residuals and their squares in
lag 1, CCF (1), is 0.04 indicating no significant correlation
and the AIC criterion favours Model II.
These results indicate that the ARFIMA(9,d,0)-
EGARCH(1,1) (Model II) model leads to further
characterisation of the HRV data. In fact, the estimate
ξˆ1 = 0.24 indicates asymmetric response to positive
and negative shocks. The conditional standard deviation
estimate σˆt represented in Fig. 1(g), captures very well the
heteroscedasticity in the original data, plotted in (a). Similar
results were obtained in other short HRV recordings.
TABLE I
ESTIMATES (STANDARD DEVIATION) FOR MODEL I AND MODEL II
FITTED TO DATA IN FIG. 1. Q(20) STANDS FOR MCLEOD-LI’ S TEST
STATISTIC FOR THE RESIDUALS.
Parameter Model I Model II
dˆ 0.34 0.34
Q(20) of ˆt 494.32 494.32
CCF(1) of ˆt,ˆ2t 0.35 0.35
uˆ0 0.13e−03(0.16e−04) −2.87 (0.20)
uˆ1 0.33 (0.02) 0.47 (0.05)
vˆ1 0.40 (0.05) 0.63 (0.03)
ξˆ1 —– 0.24 (0.02)
Q(20) of zˆt 28.84 41.40
CCF(1) of zˆt,zˆ2t 0.12 0.04
AIC of zˆt -8578.00 -8623.00
In the case of long recordings, such as ambulatory 24
hour HRV data (approximately 100000 beats), exhibiting
several non stationary characteristics, ARFIMA-EGARCH
modeling is combined with adaptive segmentation [6]: long
records are decomposed into short records of variable length
and the break points, which mark the end of consecutive
short records, are identified by AIC criterion. The short
records thus obtained have a minimum length 512 and are
subsequently modeled using ARFIMA-EGARCH models.
IV. RESULTS AND DISCUSSION
In this section, the above methodology is applied to long-
term HRV series of subjects from the Noltisalis database [9].
The results are first illustrated for the healthy subject N6 and
the patient C10, in Figs. 2 and 3, respectively.
R
R
(se
c)
0.5
1
1.5
(a)
d
0
0.5
1
(b)
Q
(2
0)
0
200
400
600
(c)
C
C
F(
1)
0
0.2
0.4
(d)
u
1,
 
v
1(-
-)
0
0.5
1
(e)
Recording time (hh:mm)
10:00 15:00 20:00 01:00 06:00 11:00
ξ
1
0
0.2
0.4
0.6
(f)
Fig. 2. (a) Tachogram of healthy subject N6, 24 hour recordings Noltisalis
database, Evolution over 24 hours of dˆ in (b), uˆ1(-) and vˆ1(- -) in (e) and
ξˆ1 in (f) estimated using ARFIMA-EGARCH models and segmentation;
statistic Q(20) of the McLeod-Li test for conditional heteroscedasticity in
(c) and cross-correlation between the residuals t and their squares of lag
1 in (d).
The long memory estimates dˆ in Fig. 2(b) and Fig. 3(b),
change over time showing circadian variation, with lowest
values during the night periods, 0 < dˆ < 0.5 in contrast
with 0.5 < dˆ < 1 for the day period. These results are in
concordance with Leite et al.[6].
The residuals from the ARFIMA modeling are tested for
conditional heteroscedasticity with McLeod-Li test. The cor-
responding test statistics Q(20) are represented in plots (c)
in Figs. 2 and 3. It is found that the percentage of segments
with heteroscedasticity are 98.8% and 94.5 % for subjects
N6 and C10, respectively. Moreover, the volatility parameters
estimates uˆ1 and vˆ1, Fig. 2(e) and Fig. 3(e), change over time
with some circadian variation for the healthy subject. For
the sick subject, the estimate vˆ1 increases and uˆ1 decreases.
Note that the estimated values for parameter v1 are over 0.5
indicating some persistence in conditional variance. These
results are in agreement with Leite et al.[6].
The residuals from the ARFIMA models are also checked
for the asymmetric effect (leverage). Plot (d) in Figs. 2 and
3 represent the values of CCF(1) between the residuals and
their squares. The correlations indicate that the conditional
variance depends also on the lagged shocks, leverage effect.
In fact, the percentage of segments with leverage effect are
93.8% and 43.8 % for subjects N6 and C10, respectively.
R
R
(se
c)
0.5
1
1.5
(a)
d
0
0.5
1 (b)
Q
(2
0)
0
200
400
600 (c)
C
C
F(
1)
0
0.2
0.4
(d)
u
1,
 
v
1(-
-)
0
0.5
1
(e)
Recording time (hh:mm)
10:00 15:00 20:00 01:00 06:00 11:00
ξ
1
0
0.2
0.4
(f)
Fig. 3. (a) Tachogram of patient C10 affected by congestive heart failure, 24
hour recordings Noltisalis database. Evolution over 24 hours of dˆ in (b), uˆ1(-
) and vˆ1(- -) in (e) and ξˆ1 in (f) estimated using ARFIMA-EGARCH models
and segmentation; statistic Q(20) of the McLeod-Li test for conditional
heteroscedasticity in (c) and cross-correlation between the residuals t and
their squares of lag 1 in (d).
Moreover, the leverage parameter estimates ξˆ1, Fig. 2(f) and
Fig. 3(f), change over time and present higher values for the
healthy subject. The same analysis is performed for patient
T3. The results for subjects N6, C10 and T3 are summarized
in Fig. 4(a).
Finally, the mean estimates dˆ, uˆ1, vˆ1 and ξˆ1 for each of
the 30 subjects in the database are summarized in Fig. 4(b).
These overall results indicate that the long memory param-
eter d is lower for healthy subjects N while the volatility
and the leverage parameters u1, ξ1 respectively, are lower
for sick subjects C and T. The parameter v1 presents high
variability for the sick subjects T. The results are promising
in differentiating health and disfunction situations deserving
further study.
V. CONCLUSIONS
This paper contributes for further characterization of the
complex dynamics of HRV with a leverage parameter. In
particular, the model under discussion indicates that values
of HRV under the mean lead to less variability in HRV than
values over the mean- asymmetric effect captured by the
leverage parameter. Furthermore this study indicates that this
effect is stronger in healthy subjects. However, further studies
are necessary to assess the importance of this parameter in
N6 C10 T3 N6 C10 T3 N6 C10 T3 N6 C10 T3
E
st
im
at
ed
-0.2
0
0.2
0.4
0.6
0.8
1
(a)
N C T N C T N C T N C T
E
st
im
at
ed
-0.2
0
0.2
0.4
0.6
0.8
1
(b)
d u
1
v
1
ξ
1
ξ
1
v
1
u
1d
Fig. 4. Boxplots for dˆ, uˆ1, vˆ1 and ξˆ1 using ARFIMA-EGARCH models
and segmentation: (a) an healthy subject-N6, a patient suffering from
congestive heart failure-C10 and a heart transplanted patient-T3, (b) subjects
mean values for each group N, C, T.
risk stratification and autonomic nervous system dysfunction
characterization.
REFERENCES
[1] Task Force of the European Society of Cardiology and North American
Society of Pacing Electrophysiology, Heart rate variability: standards
of measurement, physiological interpretation and clinical use, Circu-
lation, vol. 93, pp. 1043-1065, 1996.
[2] L. T. Mainardi, On the quantification of heart rate variability spectral
parameters using time-frequency and time-varying methods, Phil.
Trans. R. Soc. A, vol. 367, pp. 2555-275, 2009.
[3] R. T. Baillie, A. A. Cecen, and C. Erkal, Normal heartbeat series are
nonchaotic, nonlinear, and multifractal: New evidence from semipara-
metric and parametric tests, Chaos, vol. 19, pp. 028503, 2009.
[4] A. Voss, S. Schulz, R. Schroeder, M. Baumert, and P. Caminal,
Methods derived from nonlinear dynamics for analysing heart rate
variability. Phil. Trans. R. Soc. A, vol. 367, pp. 277-296, 2009.
[5] A. Leite, A. P. Rocha, and M. E. Silva, Long memory and volatility in
HRV: an ARFIMA-GARCH approach, in Computing in Cardiology,
vol. 36, pp. 165-168, 2009.
[6] A. Leite, A. P. Rocha, and M. E. Silva, Beyond long memory in
heart rate variability: an approach based on fractionally integrated
autoregressive moving average time series models with conditional
heteroscedasticity, Chaos, vol. 23, pp. 023103, 2013.
[7] M. Y. Hu, and C. Tsoukalas, Conditional volatility properties of sleep
disordered breathing, Comput. Biol. Med., vol. 36, pp. 303-312, 2006.
[8] R. S. Tsay, Analysis of Financial Time Series (2nd ed.). Wiley-
Interscience, 2005.
[9] M. G. Signorini, R. Sassi, S. Cerutti, Working on the NOLTISALIS
Database: Measurement of nonlinear properties in heart rate variability
signals. In: Proceedings of IEEE-EMBS International Conference,
Istanbul, Turkey (IEEE, Piscataway), pp. 547-550, 2001.
[10] P. Boukhan, G. Oppenheim, and M. S. Taqqu, Theory and applications
of long-range dependence, Boston: Birkha¨user, 2003.
[11] MATLAB, version 8.5 (R2015a), The MathWorks Inc., Natick, Mas-
sachusetts, 2015.


English

https://repositorio-aberto.up.pt/bitstream/10216/86116/2/157744.pdf

Modeling volatility in heat rate variability

Abstract

Similar works

Full text

Available Versions

Repositório Aberto da Universidade do Porto

Crossref

Open Repository of the University of Porto