Electronic – LTSpice AC analysis and DC analysis don’t agree

If you have a fixed frequency then you can use definite integration over one cycle to eliminate the DC. It will take one cycle for the circuit to respond, true, but you can skip that (unless you;re interested in transients, too). To avoid repeating, here's my answer giving the solution to this. I'd advise using the G+C variant, rather than the behavioural source, as the delay given from the tline is much more reliable, but the choice is yours.

If your frequency varies, then you can simply use a lowpass filter of your choice (in LTspice, avoid Laplace in .TRAN analysis). For a simple example, a 4th order Butterworth:

...or a 4th order Bessel, for almost linear phase, normalized to -3dB (see the 2.114 in .param f=1/2.114 is the frequency scaling factor):

If you DC point also varies with frequency, you're better off using the non-normalized Bessel, to have a better group delay in the passband, and to use a tline to delay your reference signal to match the Bessel delay, something like this:

Note that the way the input is delayed is just an exemplification, the way I did it, since I don't know what you have in your schematic. At the very worse, you could add a G source followed by the tline and the terminating resistor as seen in the example above.

If my way of doing is makes sense, to obtain the Laplace transforms of V(in) and V(out) which interval should I use? As you see the plots start from 0sec to 1sec. And the step ipout starts at 0.1sec. Can I just take the Laplace transform's of Vin and Vout for the entire 1 seconds and do the math? Is there any restriction here?

A better way to do this in lt spice is to do a .AC analysis, and then plot V(Vout)/V(Vin) (or whatever net names you have). AC analysis runs a sine sweep through the circuit, choose the source and the source's amplitude (the source below is set to 1V and then you can monitor any output. The green line in the below plot is \$ H(s)\$

That being said, if you really want to export time series data use the tfest command in matlab OR if you want to do system identification manually (which is not trivial and is too broad to explain all of it in an answer here, since it is a science unto itself) you could probably do this by taking the FFT of both signals to get the data into the frequency domain and then dividing the magnitude of the output by the input, you should get something that looks like a bode plot. Then use polynomial fitting to fit the data to a two pole transfer function model. (I've only done that once, usually its easier to estimate the poles by 'eyeballing' it, and much easier to use LT spice for all the computation so I only use that)

Here is the general idea (there is some translation you'll have to do from sampled to continuous and freq to laplace):

$$ \frac{|FFT(Vout)|}{|FFT(Vin)|} = |H(s)| = \frac{K}{(\tau_1 s+1)(\tau_2 s+1)}$$

$$ \frac{\angle FFT(Vout)}{\angle FFT(Vin)} = \angle H(s) = \exp^{-as}$$