New Software Kills the Burn-In Problem in Multivariate Time Series Simulation

Every statistician who has ever simulated a time series knows the ritual: generate the data, throw away the first few hundred points, and hope the remainder behaves as if it came from the true underlying process. That discarded warm-up segment, known as burn-in, has been an unavoidable tax on simulation for decades. Now a new open-source software package called Varmapack, described in the journal SoftwareX by Kristján Jónasson of the University of Iceland, eliminates the tax entirely for a broad and important class of models, and does so at speeds that dwarf the tools researchers currently rely on.

The models in question are vector autoregressive moving-average, or VARMA, processes. These extend the familiar scalar ARMA models to multiple interrelated variables at once, capturing not only how each variable depends on its own past but also how it responds to the pasts of its companions and to shared random shocks. A VARMA model of order p and q expresses the current state of the system as a weighted sum of previous states, plus a fresh Gaussian innovation, plus weighted sums of recent innovations. Such models date back to the late twentieth century and remain a workhorse of econometrics, environmental statistics, and signal analysis, often representing the same underlying dynamics more parsimoniously than pure vector autoregressive models.

The burn-in problem arises because the naive simulation recipe starts from arbitrary or zero initial values and then runs the model forward, injecting random shocks at each step. Early in the run, the series still remembers its artificial starting point, so its distribution does not yet match the true stationary distribution of the model. The standard fix is to discard an initial segment, but how long that segment must be depends delicately on the model. Highly persistent processes, whose autoregressive spectral radius approaches one, can require very long warm-ups. Jónasson’s benchmarks show the practical cost: with a default burn-in of 200 observations, the R package MTS produces a variance about seven percent too low for the first retained value of a moderately persistent test model, and comparable discrepancies of around six percent appear in Python’s Statsmodels VAR routine and Matlab’s varm when they discard the same amount.

The alternative is exact initialization, and its mathematical foundations are old. An exact algorithm for scalar ARMA simulation was published in 1978, and state-space formulations for exact VARMA simulation appeared in 1987 and 1988. The trick is to exploit the joint Gaussian distribution of the initial innovations and the initial states of the series. For a stationary model, the lagged state covariances and the mixed covariances between states and innovations can be computed in closed form, either by solving the vector Yule-Walker equations or by solving a discrete-time Lyapunov equation for the covariance of an augmented state vector. Stacking the first h states and h innovations, where h is the maximum of the autoregressive and moving-average orders, yields a joint normal distribution from which both blocks can be drawn exactly, so the very first generated term already follows the correct distribution.

Varmapack implements this idea with two complementary initialization modes. In the first, both the initial shocks and the initial states are drawn from their exact joint stationary distribution. In the second, the user supplies starting values for the series, and the package draws the corresponding initial innovations from their conditional distribution given those states, which requires the conditional covariance matrix to be positive definite. After initialization, both modes simply run the recursion forward, drawing fresh innovations and computing new states. For nonstationary models, where no stationary distribution exists, starting values must be supplied, and the package still handles the moving-average component exactly by drawing the startup innovations from their theoretical distribution subject to the constraints the model imposes.

Several technical subtleties make the implementation more than a routine coding exercise. The conditional covariance matrix can be singular even for perfectly well-behaved models; Jónasson gives the example of a two-dimensional VAR(1) model with small autoregressive coefficients and a scaled identity innovation covariance, for which the conditional covariance is singular with all entries equal to one twenty-fourth. Varmapack survives such cases by falling back on a pivoted Cholesky factorization when the ordinary factorization fails, a robustness feature borrowed from the Randompack library that supplies its random number generation. The package also corrects a small error in the earlier 2008 literature, where the recursion for mixed covariances used the wrong coefficient matrix in its final term, and it handles singular innovation covariance matrices through inverse-free formulas that remain valid in degenerate cases.

Efficiency comes from a hybrid strategy for computing the covariances that drive initialization. The vector Yule-Walker approach costs on the order of p cubed times r to the sixth power in time, where p is the autoregressive order and r the number of series, while the state-space approach, solved with the SLICOT routine SB03MD, scales with the cube of the total model order times r cubed. Neither dominates universally: the Yule-Walker route wins for small numbers of series and the state-space route for large ones, with a crossover that shifts with the model orders. Varmapack selects between them automatically using empirically determined, cross-platform cutoffs, sparing users a decision most would rather not make. The core library is written in C with level-3 BLAS operations for the heavy matrix work, and it exposes interfaces to Python, R, and Matlab, with the Python and R versions offering object-oriented model classes and the Matlab version a functional style.

The performance numbers are striking. Benchmarked on an M4 Mac with median timings per simulated value, Varmapack ran roughly 600 to 1800 times faster than MTS, six to seven times faster than the exact univariate package ts.extend, fifteen to one hundred times faster than Statsmodels VARMAX, seven to twenty-one times faster than Statsmodels VAR, and 120 to 370 times faster than Matlab’s varm. A five-platform comparison against VARMAX confirmed speed-ups of fifteen to one hundred eight times across models and architectures. Jónasson dissects these gains with dedicated analysis scripts: for packages that use burn-in, skipping it accounts for about a factor of three, with the rest split among amortizing setup costs over a thousand replicates, generating replicates in bulk, and the raw advantage of C over R or interpreted code. For VARMAX, which already avoids burn-in, the gains come roughly equally from setup amortization, bulk generation, and implementation. In one large-scale test with three hundred series, setup and generation together cost under ten nanoseconds per simulated value.

The package goes beyond simulation, computing the derived quantities researchers routinely need: theoretical autocovariances and autocorrelations, sample covariances with maximum-likelihood normalization, autoregressive and moving-average spectral radii as stationarity and invertibility diagnostics, and both standard and orthogonalized impulse response matrices obtained via Cholesky factorization of the innovation covariance. It also supports VARMAX models with exogenous inputs, which can produce nonzero, time-dependent conditional means, and simulations with fixed or time-varying mean paths. Verification is unusually thorough: test suites for all four interfaces cover twenty-three VARMA models, seeded reproducibility, error handling, and cross-checks of the C results against an independent Matlab reference implementation, including comparisons of the two covariance solvers against each other. The software is released under the MIT license, available on GitHub and through PyPI and CRAN, with planned extensions covering likelihood evaluation, estimation, structured specifications, and missing observations.

Why should anyone outside the time series community care? Because simulation with known temporal and cross-variable dependence underpins far more than academic exercises. It is central to evaluating forecasting methods, bootstrap procedures, power analyses, and teaching, and to generating long synthetic environmental and ocean-wave records from limited observations. Increasingly, it also matters for machine learning: recent deep-learning forecasters for multivariate air quality exploit exactly the kind of cross-variable structure that VARMA models encode, and controlled synthetic benchmarks generated by Varmapack could serve anomaly and fault detection research. Gains of several orders of magnitude in simulation speed translate directly into more replicates, broader parameter exploration, and larger, more reliable studies. And by removing the guesswork of choosing a burn-in length, the package removes a subtle source of bias that has quietly contaminated simulation-based comparisons for as long as researchers have been discarding their first few hundred data points.

Subject of Research: Burn-in-free exact simulation of VARMA time series with high-performance software

Article Title: Burn-in-free simulation of VARMA time series

Article References: Jónasson, K. (2026). Burn-in-free simulation of VARMA time series. SoftwareX, 36, Article 103100. https://doi.org/10.1016/j.softx.2026.103100

Image Credits: AI Generated

DOI: 10.1016/j.softx.2026.103100

Keywords: VARMA, time series, simulation, burn-in, stationary distribution, Yule-Walker equations, state-space models, open-source software, impulse response, spectral radius, Randompack, BLAS