I not sure exactly, but perhaps keep in mind definition of Vt (eg. section 3.1), Vt = Vt0 + f(Vsb), a function of source-bulk voltage. So VT03 refers to threshold voltage of M3 which has source-bulk connected (so Vt == Vt0), whereas M1 hasn't got bulk-source connected and so Vt1 depends on Vin.
Yes, large signal Ids equation and Vds>=(Vgs-Vt) is the right approach.
To consider output range - see what limits output voltage. With two-stage op-amp second stage can swing to both positive and negative rails, but distortion results if you get within one Vdsat of Vdd or Vss and M6/7 go out of saturation (so gain of the overall op-amp decreases), it depends on the application whether this is a problem. Maybe swing on output of first stage should also be considered, although as the swing here is the output swing divided by the gain of the second stage, swing here will be quite small.