> 3dREMLfit is the ideal option because fMRI time series are temporally autocorrelated
The reason is not that FMRI time series are temporally autocorrelated. Rather, we currently cannot explain everything embedded in the FMRI data, thus the resulting residuals still contain some extent of temporal structure.
> I use the -iresp command in 3dDeconvolve to model my HDR
No, the option -iresp has nothing to do the model. Instead, the model is specified through options such as stim_times, stim_files, etc. In addition, the option -iresp simply collects the results (e.g. those betas associated with each condition/task) and saves them into a separate files. Indeed, 3dREMLfit does not inherit the option -iresp from 3dDeconvolve. However, it does not matter: those betas (and their t-statistics) are directly available in the 3dREMLfit ouput and they are just buried together with other information in the overall output. So, you can still access to them.
Gang