This paper proposes a procedure for fitting Markov Modulated Poisson Processes (MMPPs)  to traffic traces that matches both the autocovariance 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 L 2-MMPPs and one M-MMPP. The 2-MMPPs are designed to match the autocovariance and the M-MMPP to match the marginal distribution. Each 2-MMPP models a specific time-scale of the data. The procedure starts by approximating the autocovariance by a weighted sum of exponential functions that model the autocovariance of the 2-MMPPs. The autocovariance tail can be adjusted to capture the long-range dependence characteristics of the traffic, up to the time-scales of interest to the system under study. The procedure then fits the M-MMPP parameters in order to match the marginal distribution, within the constraints imposed by the autocovariance matching. The number of states is also determined as part of this step. The final MMPP with M2^L states is obtained by superposing the L 2-MMPPs and the M-MMPP. We apply the inference procedure to 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, for the traces used, which include the well-known Bellcore traces.