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.