Objective: The present study proposes a model-based, statistical approach to characterizing episode patterns in paroxysmal atrial fibrillation (AF). Thanks to the rapid advancement of noninvasive monitoring technology, the proposed approach should become increasingly relevant in clinical practice. Methods: History-dependent point process modeling is employed to characterize AF episode patterns, using a novel alternating, bivariate Hawkes self-exciting model. In addition, a modified version of a recently proposed statistical model to simulate AF progression throughout a lifetime is considered, involving non-Markovian rhythm switching and survival functions. For each model, the maximum likelihood estimator is derived and used to find the model parameters from observed data. Results: Using three databases with a total of 59 long-term ECG recordings, the goodness-of-fit analysis demonstrates that the proposed alternating, bivariate Hawkes model fits SR-to-AF transitions in 40 recordings and AF-to-SR transitions in 51; the corresponding numbers for the AF model with non-Markovian rhythm switching are 40 and 11, respectively. Moreover, the results indicate that the model parameters related to AF episode clustering, i.e., aggregation of temporal AF episodes, provide information complementary to the well-known clinical parameter AF burden. Conclusion: Point process modeling provides a detailed characterization of the occurrence pattern of AF episodes that may improve the understanding of arrhythmia progression.