Extracting phase-shift and gain from a time series information
This seems like a great chance to use the discrete fourier series.  It is possible to avoid all the nasty transformations and only pull one term from the series, and correlate your data with only one desired frequency for comparison.  This method essentially uses only one term of the discrete fourier series to correlate the timeseries data with only one discretized frequency.  The math gets pretty nasty, and I'm not great with the mathscript on this site yet, so I'll provide you an external link to an eloquent explanation of the theory behind it.  I'll warn you, it's very involved, and frankly unnecessary for simple applications.

https://www.math.ubc.ca/~feldman/m267/dft.pdf

I myself will provide some application examples with matlab code, easily ported to something like python if you don't have access.

What we need to do is take the time series data, and multiply each data point by the value of a cosine wave at that time point, and sum them all up.  Then repeat with a sin wave.  Yep, all that nasty math, and it boils down to something that simple.  That will give us coefficients for the correlation between our time series and the sin and cosine waves of the desired frequency.  Once we have these two numbers, use the standard magnitude and phase equations:

mag = sqrt(a^2+b^2)

phase = atan2(a,b)

Here is some code that generates a sine wave, simulates an output from an arbitrary linear system, and compares their magnitude and phase.  I set up my data to have 1000 data points and a timestep of .001, you will need to change these parameters to match your sample rates:

N=1000; %simulation parameters, adjust to your measurement frequency
T = 0.001;

freq=5*2*pi; %5 Hz

n = 0:1:(N-1);
t = n*T;


y=sin(freq*t);

G = tf(100,[1 100]); %arbitrary system for example

y1=lsim(G,y,t)'; %get a linear simulation output from the example system

plot(t,y1,t,y)



a_sum=0;
b_sum=0;
w = 5*2*pi; %declare relevant frequency for comparison
for (k=0:(N-1)) %iterate over all samples
    time = k*T;
    a_sum = a_sum + y(k+1)*cos(w*time);%perform fourier correlation with ONLY THE RELEVANT FREQUENCY
    b_sum = b_sum + y(k+1)*sin(w*time);
end
a = a_sum*2/N;
b = b_sum*2/N;
mag = sqrt(a^2+b^2)  %magnitude formula
phase = atan2(a,b)  %phase is arctangent of fourier coefficients


a_sum=0;  %now do it all again for the output wave
b_sum=0;
w = 5*2*pi;
for (k=0:(N-1))
    time = k*T;
    a_sum = a_sum + y1(k+1)*cos(w*time);
    b_sum = b_sum + y1(k+1)*sin(w*time);
end
a = a_sum*2/N;
b = b_sum*2/N;
mag1 = sqrt(a^2+b^2)
phase1 = atan2(a,b)


code outputs:



mag =

    1.0000


phase =

  -6.5239e-17


mag1 =

    0.9539


phase1 =

   -0.2984


As you can see, I generated a 5Hz (or 5*2*pi radian) sine wave, simulated an output not unlike the one you see in your time series, and performed the necessary correlations in for loops to compute the magnitude and phase.  To get the magnitude and phase shift you only need to take the difference between the two magnitude and phase outputs, and boom, you've got it.

Let me know if I can edit my answer to clear anything up, I've crammed a lot of discrete-time systems analysis material into a very small space here so I'm happy to expand if necessary.