I don't see the motive for doing this. There are 2 time series regression programs involved.
3dDeconvolve does regression using ordinary least squares (OLSQ).
3dREMLfit does regression using generalized least squares (GLSQ), with REML estimation of the covariance matrix involved.
Your script would do GLSQ estimation of the motion effects and subtract those from the data. Then with that subtracted data, you would then do OLSQ estimation of the task parameters.
I simply can't see why. It is certainly
possible, but I don't think it is
reasonable. If you believe the noise in the time series is correlated, then why would you seek statistical measures of the task parameters using OLSQ and of the motion parameters using GLSQ?