Thank you Rick for your input!
Like you have pointed out, I have been using 3dmaskdump and 3dUndump for this purpose.
I understand your concern of high correlation values. In my analysis, I use a spatial mask of 19 voxels over a "3d+t" of 463 time points, which I think is OK. The issue I have here is the computational load since I need to move the mask in the brain as a "search light", so the procedure of "3dmaskdump + Least-squre-fit" has to be done for each voxel; and what's even crazier is to build a null distribution of this correlation by shuffling the original time series for like 1000 times. I just wish that I had options of "-canocorr" and "-jobs" in "3dTcorrelate". ^_^

Thank you again for your help!