This field of control systems is called System Identification. If you have access to it, the System Identification Toolbox from Matlab is precisely what you are looking for. The associated documentation should also tell you what you need to know about this subject. This is also a common topic in Control Systems Textbooks. Skip through any you have available and see if any have a section on System Identification.

Dear Varun,

what you want to do is essentially to fit a rational complex function of a chosen order to a set of complex data (frequency response points). A natural requirement is to obtain a fitting function with minimal distance to your data points at the given frequencies. This leads to a nonlinear least squares optimization problem which has to be solved numerically. The Levenberg-Marquardt method is an example of a robust algorithm which can be used for this purpose (see e.g. the publication The Levenberg-Marquardt method for nonlinear least squares curve-fitting problems by H. P. Gavin).

The given problem can be approximated by an autoregressive model leading to classical linear least squares (see the paper A poly-reference implementation of the least-squares complex frequency-domain estimator by Guillaume et. al). It will generally produce biased estimates of the model parameters. However, they may serve as an initial condition for the second step of nonlinear optimization. The results will be close to the optimum in the case of negligible noise in your frequency response data. Another interesting approach is the method of orthogonal least squares which is also able to estimate the necessary order of the model transfer function (Billings and Tsang, An orthogonal estimation algorithm for complex number systems).

Numerical robustness of the least squares estimators is an issue due to the ill conditioning of the problem. Selection of a proper polynomial basis functions may be needed to get suitable results (see e.g. the paper On numerically reliable frequency-domain system identification by Voorhoeve et al.)

The above mentioned System identification toolbox in Matlab provides a set of suitable routines.

The area which would address your doubt would be the area of System identification & Parameter estimation broadly. But you would be more concerned about the system identification part . I would suggest to go through the text book " system identification theory for the user" by Lennart Ljung for the basic methods of system identification. But if you want to have a quick and real time solution for your doubt you can go with MATLAB system identification toolbox also. 

transfer function can be simply obtained by identifying the corner frequncies and low frequency magnitude point and refer I J Nagrath and gopal

And also examine the phase plot changes  at corner frequencies

Keep things simple and assume a 2nd-order plant. If it looks like it is a low-pass plant, then place a pair of repeated poles on the (-ve) real axis. To increase the low-pass character, slide the poles along the real axis towards imaginary axis. If it looks like your plant has resonant peaks away from dc, then use a complex conjugate pair of poles. Set the imaginary coordinate of the poles, in the left-hand side of the s-plane, to match the resonant frequency; slide the poles towards the imaginary axis, to increase the height of the peak. 

If the response is an asymptotic plot then we can easily obtain by basic corner frequency approach.But here we have data points only and that too its not an asymptotic plot. So the basic corner frequency approach may not come handy here.

It depends on how noisy or clean your data is. If it is noisy, then you may need to apply some filter to it based on known information about the data. Examples could be Kalman filtering, Baisian filters, or high, low pass filters of 2nd, 3rd or higher orders. Once the data is clean, you need to answer the question: do I know the input or the output, both, none, or do I know about them. based on each answer there is different technique to solve the problem. there are some new techniques that use wavelets, see my paper, Wavelets as a tool for systems analysis and control. Usually the result you get is in frequency domain (on the jw axis) to convert to s domain, you need to substitute every jw with s, and that would be the transfer function. 

Dear Varun,

E. O. Doebelin 

https://www.amazon.com/System-Modeling-Response-Ernest-Doebelin/dp/0471032115/ref=sr_1_1?ie=UTF8&qid=1470685278&sr=8-1&keywords=System+modeling+and+response

describes Levy´s method. This method takes frequency response data (real and imaginary values at some frequencies) and returns the parameters of a transfer function in the s domain. The user must choose the model structure. The procedure is based on "approximage Least Squares" and is very easy to implement. The only recommendation is not to use infomation much above the system pass-band. The reason for this (I assume) is that the phase at medium-high frequencies is very inaccurate. Apart from that, the method works well.