Cerebrovascular reactivity (CVR) can be mapped noninvasively using blood oxygenation level dependent (BOLD) fMRI during a breath-hold (BH) task. Previous studies showed that the BH BOLD response is best modeled as the convolution of the partial pressure of end-tidal CO2 (PetCO2) with a canonical hemodynamic response function (HRF). However, previous model comparisons employed a global bulk time lag, which is now well accepted to provide only a rough approximation of the heterogeneous distribution of response latencies across the brain. Here, we investigate the best modeling approach for mapping CVR based on BH BOLD-fMRI data, when using a lagged general linear model approach for voxelwise lag optimization. In a group of 14 healthy participants, we compared two types of regressors (PetCO2 and Block), and three convolution models (no convolution; convolution with a single gamma HRF; and convolution with a double gamma HRF), as well as a range of HRF delays and dispersions (for models with convolution). Convolution with a single gamma HRF yielded the greatest CVR values in PetCO2 models, while a double gamma HRF performed better for block models. Although PetCO2-based regressors generally outperformed block-based regressors, as expected, the latter may be an appropriate alternative in cases of poor CO recordings. Overall, our results support the use of specific modeling approaches for CVR mapping based on end-expiration BH BOLD-fMRI, including the voxelwise optimization of the lag.