Pulses from radio pulsars become dispersed as they pass through the ionized interstellar medium. This dispersive effect delays pulses as a function of radio frequency as well as the DM, the integrated line-of-sight electron density. Changes in the Earth-pulsar line of sight cause well-known, significant time variability of the DM. We have developed the statistical framework to describe DM timeseries as Seasonal (S) Continuous (C) Autoregressive (AR) Fractionally (F) Integrated (I) Moving Average (MA) processes. The advantage of the SCARFIMA approach is that it can naturally separate the various contributions to the DM variations and that since we understand the underlying physics, we do not have to identify the order of the ARMA component via complex model selection. The turbulent component of the medium produces variations described by a power-law red-noise spectrum whose index can be estimated with an ARFIMA(0, d, 0) process after differencing. This differencing also naturally removes contamination of the timeseries by any linear trends in the DM due to radial Earth-pulsar motions, which can then be recovered Priors on the amplitude of the spectrum can be made using pulsar scintillation measurements. Variations due to the line of sight crossing the solar wind can be modeled as a SARMA process with a yearly timescale and low-order ARMA components. The continuous analogue to SARFIMA is required for the irregularly sampled timeseries with variable measurement errors. The maximum likelihood approach to parameter estimates can be extended to multiple DM timeseries to create a pulsar interstellar medium array (PISMA) with the goal of measuring correlated signatures in pulsar data due to the interstellar medium and solar wind, analogously to the observation of a pulsar timing array (PTA) with the goal of measuring correlated signatures due to gravitational waves.