Hi,
Consider the following situation:

      Take a MOS transistor and give  Vdc+ v.sin(wt) to the source(S)  and Vdc-v.sin(wt) to the drain(D).
Fix the gate(G)  potential at Vg, so that the MOS transistor is biased in the linear region, i.e, Vgs > Vth and Vds=0 for a DC operating point. Do a tran analysis for a sufficiently long time and measure the third harmonic current at the drain(Id3) by FFT analysis.
        The observation is that the magnitude of Id3 just grows as the square of v, but we expect it to grow as cube of v.

For my case, I'm using UMC 0.13u trasnsistors, with Vdc=150 mV, Vg=600 mV , v=1 mV and frequency= 3 MHz.
I tried both with Spectre (version: 5.0.0) and Hspice (version: W-2005.03-SP1 ), but I'm getting the same results.
 
 Can someone explain why Id3 is not proportional to v cubed, as expected?

The Ids(Vds) characteristic around Vds=0 is not modeled very well in BSIM3. There is nothing in the model to make sure that the derivative of the Ids(Vds) characteristic is continuous (and thus equal to zero) at Vds=0. So, the Ids(Vds) characteristic of the BSIM3 model usually has an unphysical kink at Vds=0 which causes all sorts of strange and unrealistic behavior in simulation.

Yes, you should probably try that. I have not used EKV myself so far (most foundries only support BSIM3) but as far as I know, it does not have this problem.

By the way, I just realized that my statement above was not quite correct. It is not the Ids(Vds) characteristic that has the kink but the conductance as a function of Vds. Due to symmetry, the even-order derivatives of Ids with respect to Vds must be zero at Vds=0. This is generally not the case in BSIM3 for short-channel devices. This subject is treated in some more detail in section 4.7 of the book by William Liu (see http://www.designers-guide.org/Books/#Liu-2001).