Investigating the relationship between task-related hemodynamic responses and cortical excitability is challenging because it requires simultaneous measurement of hemodynamic responses while applying noninvasive brain stimulation. Moreover, cortical excitability and task-related hemodynamic responses are both associated with inter-/intra-subject variability. To reliably assess such a relationship, we applied hierarchical Bayesian modeling. This study involved 16 healthy subjects who underwent simultaneous Paired Associative Stimulation (PAS10, PAS25, Sham) while monitoring brain activity using functional Near-Infrared Spectroscopy (fNIRS), targeting the primary motor cortex (M1). Cortical excitability was measured by Motor Evoked Potentials (MEPs), and the motor task-related hemodynamic responses were measured using fNIRS 3D reconstructions. We constructed three models to investigate: (1) PAS effects on the M1 excitability, (2) PAS effects on fNIRS hemodynamic responses to a finger tapping task, and (3) the correlation between PAS effects on M1 excitability and PAS effects on task-related hemodynamic responses. Significant increase in cortical excitability was found following PAS25, whereas a small reduction of the cortical excitability was shown after PAS10 and a subtle increase occurred after sham. Both HbO and HbR absolute amplitudes increased after PAS25 and decreased after PAS10. The probability of the positive correlation between modulation of cortical excitability and hemodynamic activity was 0.77 for HbO and 0.79 for HbR. We demonstrated that PAS stimulation modulates task-related cortical hemodynamic responses in addition to M1 excitability. Moreover, the positive correlation between