This paper proposes a fitting procedure for Markov Modulated Poisson Processes (MMPPs) that matches both the autocovariance tail and marginal distribution of the counting process. A major feature of the procedure is that the number of states is not fixed a priori. It is an output of the fitting process thus allowing the number of states to be adapted to the particular trace being modeled. The MMPP is constructed as a superposition of one M/2-MMPP and one 2-MMPP. The M/2-MMPP is designed to match the density function and the 2-MMPP to match the autocovariance tail. The procedure starts by approximating the autocovariance by a weighted sum of two exponential functions. The exponential with slower decay is selected to model the autocovariance of the 2-MMPP. The autocovariance tail can be adjusted to capture the long-range dependence characteristics of the traffic. The procedure then fits the M/2-MMPP parameters in order to match the density function, within the constraints imposed by the autocovariance matching. The number of states is also determined as part of this step. The final MMPP with M states is obtained by superposing the 2-MMPP and the M/2-MMPP. We apply the inference procedure to IP traffic traces exhibiting long-range dependence and evaluate its queuing behavior through simulation. Very good results are obtained, both in terms of queuing behavior and number of states, in particular, with the well-known Bellcore traces.