There are 3 sources of estimation error:

(a) the hemodynamic response function (HRF or IRF) model is inappropriate -- not much we can do about that, so we just ignore it

(b) intra-subject noise; the standard deviation for the HRF from this is output using the -sresp option in 3dDeconvolve

(c) inter-subject variability; this effect is probably larger than (b), unless you had very few task repetitions per subject

If we ignore (b), then you could use 3dcalc to calculate the standard deviation of each voxel at each lag, and divide that value by the square root of the number of subjects; the result would be the SEM of the HRF, which you could then use to plot error bands about the grand mean HRF.