Modelling the amplitude and phase through an amplifier
I apologize, the following is a bit rough, but I've put together an example in MATLAB illustrating all the different elements I think you need. That is: 


a fist-order model for an amplifier (a common approximation for most op. amps.)
how to discretize it to find an IIR filter
one way to implement a realization of the IIR filter (straight forward state-space implementation)
how to generate a transfer-function estimate from measured data
one way (not the most robust, but perhaps the most obvious) to fit the transfer-function coefficients to the transfer-function estimate


To read up on the theory behind all this I would recommend the following text-books:

Chi-Tsong Chen, Linear System Theory and Design (3rd edition). Oxford University Press, 1999.

Petre Stoica and Randolph Moses, Spectral Analysis of Signals. Prentice Hall, 2004.

Lennart Ljung, System Identification: Theory for the User (2nd edition). Prentice Hall, 1999.

The last book is somewhat heavy, but concentrate on "output error" methods, and you should be very well covered for systems where you know the order of the model, or "subspace" methods if you don't. These methods are available in the MATLAB system identification toolbox.

%% Amplifier model
wb = 1e4; % Amplifier bandwidth
ga = 1; % Amplifier gain

% Modeling amplifier as having fist-order dynamics
% If using Octave run these in the command window:
%pkg install -forge control
%pkg load control
H = tf(ga,[1/wb 1]); % Continuous-time transfer function

%% Finding a realization (convert to state-space formulation)
% V/U = ga/(s/wb + 1)
% V*(s/wb + 1) = ga*U
% s*V = -wb*V + ga*wb*U
% dv/dt = -wb*v + ga*wb*u
%
% Standard form:
% dx/dt = A*x + B*u
% y = C*x + D*u
%
% with x = v, y = x
A = -wb; B = ga*wb; C = 1; D = 0;

%% Finding a discretization:
Ts = 1/1e5; % sampling time
% Using Euler integration
Ad1 = 1 + Ts*A; Bd1 = Ts*B; Cd1 = C; Dd1 = D;
% Using exact discretization (use when system has analog-digital converters)
Ad2 = exp(A*Ts); Bd2 = (Ad2 - 1)/A*B; Cd2 = C; Dd2 = D;

%% Computing a step-response to demonstrate one way to implement an IIR filter
u = [zeros(1,100) ones(1,200)]; % input
x = zeros(2,numel(u)+1);
y = zeros(2,numel(u));

t = 0:Ts:(numel(u)-1)*Ts;

% Comupting filter outputs
for k = 1:numel(u)
    x(1,k+1) = Ad1*x(1,k) + Bd1*u(k); % response of Euler integration scheme
    y(1,k) = Cd1*x(1,k) + Dd1*u(k);

    x(2,k+1) = Ad1*x(2,k) + Bd1*u(k); % response of exact discretization
    y(2,k) = Cd1*x(2,k) + Dd1*u(k);
end

figure(1)
plot(t,y(1,:),t,y(2,:))

%% Computing the response to noise for identification of parameters
% The frequency response of a system is in general most reliably found when
% exciting the system with white noise. A pseudo-random binary sequence is
% theoretically optimal when you don't know anything about your system, but
% can in practice cause problems due to the switched nature. A Gaussian
% noise source is usually a good substitute.
un = randn(1,1e4); % input
xn = zeros(1,numel(un)+1);
yn = zeros(1,numel(un));

% Comupting filter outputs
for k = 1:numel(un)
    xn(k+1) = Ad2*xn(k) + Bd2*un(k);
    yn(k) = Cd2*xn(k) + Dd2*un(k);
end
% Remove transient as you want to look at stationary response
un_ = un(201:end);
yn_ = yn(201:end);

% Compute power spectral densitites and cross-power spectral densitites
% using a simplified version of the Welch method
N = length(yn_); % length of original sequence
L = 5; % number of PSD averages
M = floor(N/L); % length of sequence segments
W = 0:(2*pi)/M:2*pi-(2*pi)/M; % normalized PSD frequencies

WIN = hanning(M); WIN = WIN(:).'; % optional windowing (if the segments are short)
Pwin = 1/M*sum(abs(WIN).^2); % window "power" correction (for Welch method)

% PSD for input sequence
Puu = zeros(1,M);
for k = 1:L
    un_segment = un_((k-1)*M+1:k*M);
    un_segment = un_segment(:).';
    un_windowed = un_segment.*WIN;
    Uft = 1/sqrt(2*pi*M)*fft(un_windowed);

    Puu = Puu + (Uft.*conj(Uft));
end
Puu = 1/(L*Pwin)*Puu;

% Cross-PSD between input and output sequence
Puy = zeros(1,M);
for k = 1:L
    un_segment = un_((k-1)*M+1:k*M);
    un_segment = un_segment(:).';

    yn_segment = yn_((k-1)*M+1:k*M);
    yn_segment = yn_segment(:).';

    yn_windowed = yn_segment.*WIN;
    Yft = 1/sqrt(2*pi*M)*fft(yn_windowed);

    un_windowed = un_segment.*WIN;
    Uft = 1/sqrt(2*pi*M)*fft(un_windowed);

    Puy = Puy + (Yft.*conj(Uft));
end
Puy = 1/(L*Pwin)*Puy;

% Transfer-function estimate
Huy = Puy./Puu;

figure(2)
subplot(2,1,1)
semilogx(W,abs(Huy))
subplot(2,1,2)
semilogx(W,unwrap(atan2(imag(Huy),real(Huy))))

%% Fitting a model to the output data

% finding the z-transform of discrete-time model (Euler integration)
% using the standard expression
z = tf('z',Ts);
Hd1 = Cd1*(z - Ad1)\Bd1 + Dd1; 

% finding the z-transform of discrete-time model (exact discretization)
% using an alternative method
b0 = Bd2*Cd2;
a0 = -Ad2;
a1 = 1;
Hd2 = tf(b0,[a1 a0],Ts);

% compare the continuous and discrete-time models
figure(3)
bode(H,Hd1,Hd2)

zv = exp(1i*W); % starred transform (between z- and s-domains)

% H = b0/(a1*z + a0)
% H(W) = b0/(a1*exp(i*W) + a0)
% a1 = 1; % in order to have unique coefficients for a0 and b0
% H(W) = b0/(exp(i*W) + a0)
%
% (exp(i*W) + a0)*H(W) = b0
% exp(i*W)*H(W) = -a0*H(W) + b0
% exp(i*Ts*w)*H(W) = [-H(W); 1]*[a0, b0]
% zeta = phi*theta;

zeta = zv(:).*Huy(:);
phi = [-Huy(:), ones(size(Huy(:)))];

theta_true = [a0; b0]
theta_est = phi\zeta