Transfer function estimation using system identification
Let us denote the output feedback SISO plant transfer function $\frac{Y(s)}{U(s)}=P(s)$ so that the open-loop (OL) transfer function of the system loop is $CP$ where $C(s)=K_p + K_d s + \frac{K_i}{s}$ and the closed-loop (CL) transfer function is $\frac{CP}{1+CP}$. From the attempt at the solution to the problem in the OP, we can see that a second order model might be considered a good idea.
Since the controller transfer function ($K_p, K_d s, {K_i}$) and the system is using direct output feedback, the control input (or actions) may be obtained from the provided data as $u(t) = K_p y(t) + K_d \dot{y}(t) + K_i \int_0^t y(t) \; dt$. The data available is then $\{y(t_i), u(t_i)\}$ which can be directly fed into a system identification tool. Alternately one may used least squares to obtain the second order model parameters ($A, B)$ to fit the state transition model ${\Huge[} \begin{matrix} y \\ \dot{y}\end{matrix}{\Huge]} = A {\Huge[} \begin{matrix} y \\ \dot{y}\end{matrix}{\Huge]} + B {\Huge[} \begin{matrix} 0 \\ u\end{matrix}{\Huge]}$. Care must be taken not to use the CL data to fit the model instead since, as clearly seen from the model, there is no integral term in the state space representation the model is attempting to use. Finally, from a preliminary observation of the graphical representation of the data, the data can be divided into three distinct data set split by time-intervals in which the system undergoes large control actuations which seem to become saturated. These three data sets are observed to have the typical second order behavior reflecting a first or second order system controlled by a PID controller.