This paper proposes a parameter fitting procedure using Markov Modulated Poisson Processes (MMPPs) that leads to accurate estimates of queuing behavior for network traffic exhibiting long-range dependence behavior. The procedure matches both the autocovariance and marginal distribution of the counting process. A major feature is that the number of states is not fixed a priori, and can 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.