Paper deep dive
QARIMA: A Quantum Approach To Classical Time Series Analysis
Nishikanta Mohanty, Bikash K. Behera, Badshah Mukherjee, Pravat Dash
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 95%
Last extracted: 4/10/2026, 4:28:53 AM
Summary
QARIMA is a hybrid quantum-classical time series analysis framework that replaces classical ARIMA components with quantum-inspired algorithms. It utilizes swap-test-driven quantum autocorrelation (QACF) and quantum partial autocorrelation (QPACF) for lag discovery, and employs fixed-configuration Variational Quantum Circuits (VQCs) for parameter estimation (VQC-AR and VQC-MA) and weak-lag refinement, aiming to reduce meta-optimization overhead while maintaining interpretability.
Entities (6)
Relation Signals (3)
QARIMA → integrates → Variational Quantum Circuits
confidence 95% · integrates quantum-assisted lag discovery with fixed-configuration variational quantum circuits (VQCs)
QARIMA → uses → Compact Swap Test
confidence 95% · Differencing and candidate lags are identified via swap-test-driven quantum autocorrelation
QACF → identifies → candidate lags
confidence 90% · Differencing and candidate lags are identified via swap-test-driven quantum autocorrelation (QACF)
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We present a quantum-inspired ARIMA methodology that integrates quantum-assisted lag discovery with \emph{fixed-configuration} variational quantum circuits (VQCs) for parameter estimation and weak-lag refinement. Differencing and candidate lags are identified via swap-test-driven quantum autocorrelation (QACF) and quantum partial autocorrelation (QPACF), with a delayed-matrix construction that aligns quantum projections to time-domain regressors, followed by standard information-criterion parsimony. Given the screened orders $(p,d,q)$, we retain a fixed VQC ansatz, optimizer, and training budget, preventing hyperparameter leakage, and deploy the circuit in two estimation roles: VQC-AR for autoregressive coefficients and VQC-MA for moving-average coefficients. Between screening and estimation, a lightweight VQC weak-lag refinement re-weights or prunes screened AR lags without altering $(p,d,q)$. Across environmental and industrial datasets, we perform rolling-origin evaluations against automated classical ARIMA, reporting out-of-sample mean squared error (MSE), mean absolute percentage error (MAPE), and Diebold--Mariano tests on MSE and MAE. Empirically, the seven quantum contributions -- (1) differencing selection, (2) QACF, (3) QPACF, (4) swap-test primitives with delayed-matrix construction, (5) VQC-AR, (6) VQC weak-lag refinement, and (7) VQC-MA -- collectively reduce meta-optimization overhead and make explicit where quantum effects enter order discovery, lag refinement, and AR/MA parameter estimation.
Tags
Links
- Source: https://arxiv.org/abs/2604.08277v1
- Canonical: https://arxiv.org/abs/2604.08277v1
Trouble viewing inline? Open PDF directly →
Full Text
89,450 characters extracted from source content.
Expand or collapse full text
QARIMA: A Quantum Approach To Classical Time Series Analysis Nishikanta Mohanty 1* , Bikash K. Behera 2,3 , Badshah Mukherjee 4 , Pravat Dash 5 1 Centre for Quantum Software and Information, University of Technology Sydney, 15 Broadway, Ultimo, Sydney, 2007, NSW, Australia. 2 Bikash’s Quantum (OPC) Pvt. Ltd., Balindi, Mohanpur, 741246, WB, India. 3 Universit`a degli Studi di Cagliari, Via Is Mirrions, Cagliari, 09123, Italy. 4 SAS Middle East Fz-Llc., Dubai, 9262, UAE. *Corresponding author(s). E-mail(s): nishikanta.mohanty@uts.edu.au; Contributing authors: bikas.riki@gmail.com; badshah.mukherjee@outlook.com; pravat.dash@outlook.com; Abstract We present a quantum-inspired ARIMA methodology that integrates quantum- assisted lag discovery with fixed-configuration variational quantum circuits (VQCs) for parameter estimation and weak-lag refinement. Differencing and can- didate lags are identified via swap-test–driven quantum autocorrelation (QACF) and quantum partial autocorrelation (QPACF), with a delayed-matrix construc- tion that aligns quantum projections to time-domain regressors, followed by standard information-criterion parsimony. Given the screened orders (p, d, q), we retain a fixed VQC ansatz, optimizer, and training budget preventing hyper- parameter leakage and deploy the circuit in two estimation roles: VQC-AR for autoregressive coefficients and VQC-MA for moving-average coefficients. Between screening and estimation, a lightweight VQC weak-lag refinement re-weights or prunes screened AR lags without altering (p, d, q). Across environmental and industrial datasets, we perform rolling-origin evaluations against automated classical ARIMA, reporting out-of-sample mean squared error (MSE), mean absolute percentage error (MAPE), and Diebold–Mariano tests on MSE and MAE. Empirically, the seven quantum contributions (1) differencing selection, (2) QACF, (3) QPACF, (4) swap-test primitives with delayed-matrix construction, (5) VQC-AR, (6) VQC weak-lag refinement, and (7) VQC-MA collectively reduce 1 arXiv:2604.08277v1 [quant-ph] 9 Apr 2026 meta-optimization overhead and make explicit where quantum effects enter order discovery, lag refinement, and AR/MA parameter estimation. Keywords: Compact Swaptest, Variational Quantum Circuits, Quantum Optimisation, ARIMA, Quantum-ACF, Quantum-PACF, Forecast Error Metrics, Diebold-Mariano test 1 Introduction Time series data forecasting is an essential topic in statistical modelling and machine learning, with applications across economics, energy, healthcare, climate research, and more. The AutoRegressive Integrated Moving Average (ARIMA) framework [1] is one of the most prevalent models for both stationary and non-stationary univari- ate sequences, since it decomposes signals into Auto-Regressive (AR), Differencing (I), and Moving-Average (MA) components. Although ARIMA models are inter- pretable and computationally efficient, they often fall short in capturing long-range dependencies, nonlinear dynamics, or stochastic effects that are inherent in real- world situations. Furthermore, ARIMA significantly depends on traditional heuristics, including partial autocorrelation functions (PACF), information criteria (e.g., AIC, BIC), and residual-based diagnostics for order estimation and parameter optimization [2]. These approaches may become fragile or computationally intensive in high- dimensional or noisy environments. Recent advancements in quantum computing and quantum-inspired algorithms have created opportunities to overcome these limitations by reconsidering classical learning objectives via quantum principles [3]. Quantum sim- ilarity measures (e.g., swap tests, quantum cosine projections), variational quantum circuits (VQCs), and hybrid optimisation loops have shown potential in enhancing learning for structured data [4, 5]. Nonetheless, a comprehensive integration of these methods into time series modelling, especially within the ARIMA framework, remains little investigated. We introduce a hybrid, quantum–inspired ARIMA framework that replaces key classical steps with quantum counterparts while preserving ARIMA’s interpretability and modularity. The pipeline comprises seven coordinated components: (1) differ- encing selection, (2) quantum autocorrelation (Q–ACF), and (3) quantum partial autocorrelation (Q–PACF) obtained from swap–test projections, together with (4) a delayed–matrix construction aligning quantum measurements to time–domain regres- sors. Given the screened orders (p,d,q), estimation proceeds with a fixed–configuration VQC (fixed ansatz, optimizer, and budget) used in two roles: (5) VQC–AR for autoregressive coefficients b and (7) VQC–MA for moving–average coefficients θ, each optimized under a norm–constrained, quantum–augmented loss. Between screening and estimation, (6) a lightweight VQC weak–lag refinement re–weights or prunes screened AR lags without altering (p,d,q). Instead of selecting lags via classical PACF/ACF heuristics or sweeping AIC/BIC grids, we estimate QACF and QPACF using compact-swap–test projections and phase- corrected cosine similarities; candidate (p,q) are gated by these quantum diagnostics, 2 and only then do we apply information-criterion parsimony. This keeps order discovery consistent with the projection logic used later in VQC-refined AR/MA estimation and weak-lag pruning. Classical ARIMA treats model selection and estimation as largely separate, our system integrates order discovery, weak–lag refinement, and parameter estimation within a unified quantum–classical loop that maintains an interpretable ARIMA form. The resulting methodology is theoretically grounded, scalable across datasets, and implementable on near-term hardware through simulation backends, providing a principled bridge between statistical forecasting and quantum machine- learning techniques. 2 Background ARIMA models are among the most widely used approaches for modeling and forecast- ing univariate time series. The method was originally formalized by Box and Jenkins [6] and has since become foundational in econometrics and applied statistics. Given a time series y t N t=1 , an ARIMA (p,d,q) process is defined by differencing the original series d times to induce stationarity and then modeling the result as a combination of AR and MA components. Formally, the model is expressed as: B p (L) (1− L) d y t = Θ q (L)ε t ,(1) where: • L is the lag operator, Ly t = y t−1 . • B p (L) = 1− b 1 L−·− b p L p denotes the AR polynomial. • Θ q (L) = 1 + θ 1 L +· + θ q L q denotes the MA polynomial. • ε t is white noise, typically assumed Gaussian, ε t ∼N (0,σ 2 ). After applying (1− L) d differencing, the ARIMA model specifies that the current observation depends linearly on its past p values and on q past forecast errors: y t = p X i=1 b i y t−i + ε t + q X j=1 θ j ε t−j .(2) ARIMA models are commonly applied for: • Forecasting economic indicators (e.g., GDP, inflation). • Modeling financial time series (e.g., stock prices, returns). • Predicting demand in supply chain applications. Parameters (p,d,q) are typically identified via criteria such as AIC, BIC, or cross- validation. In our preceding section of 3, we will detail the process of calculating the parameters (p,d,q) using a combination of quantum processes and classical optimiza- tion. These parameters will be further tuned by a VQC for improved performance in case-to-case scenarios. 3 3 Proposed Methodology Our methodology proposes a modular approach to construct Quantum-ARIMA that consists of several distinct phases, each complementing the other. The workflow begins by preparing a lagged representation of the series and selecting the differencing order required for approximate stationarity. Autocorrelation and partial–autocorrelation structure are then inferred by quantum projections: compact swap–test primitives provide phase–corrected cosine similarities that we aggregate into quantum analogues of ACF and PACF. These diagnostics are computed from compact swap-test projec- tions with phase-corrected cosine similarities and feed directly into the delayed-matrix construction used for estimation. A delayed–matrix construction then aligns the quan- tum measurements with time–domain regressors, and standard information criteria are subsequently applied to retain a concise order triple (p ⋆ ,d ⋆ ,q ⋆ ). Given these screened orders, estimation proceeds with the same fixed VQC configuration at every rolling origin. Autoregressive coefficients φ are obtained by minimizing a quantum–augmented loss that combines a prediction–error term with cosine–alignment and projection–entropy penalties derived from the swap–test out- puts; an ℓ 2 constraint stabilizes the solution. Between screening and estimation we introduce a lightweight weak–lag refinement that re–weights or prunes the screened AR lags without altering (p ⋆ ,d ⋆ ,q ⋆ ), allowing the model to discourage marginal lags while preserving the selected order. Moving–average coefficients θ are then estimated from the delayed innovation matrix using the same loss family; initialization follows a conditional–least–squares (OLS) fit on the innovations with coefficients clipped for numerical stability. Because the circuit, optimizer, and budget are held constant, the effect of the variational layer is comparable across datasets and time. This section provides the architectural overview; the next section details the con- crete procedures. The differencing routine, quantum ACF/PACF construction, and delayed–matrix alignment appear as Algorithms 4–8. The AR estimator, weak–lag refinement, and MA estimator are specified in Algorithms 10–17, and the final ARMA composition used for forecasting is summarized in Algorithm 18. Together these components deliver an interpretable ARIMA model in which quantum effects enter precisely at order discovery, weak–lag refinement, and AR/MA estimation, with all other choices fixed for reproducibility. 4 Processes and Algorithms In this section, we will detail the algorithms that are used to calculate Quantum ARIMA components for parameter estimation. Along with these, we will describe some supporting mechanisms that are used to aid the major algorithms. 4.1 Quantum-Inspired Cosine Similarity via Compact Swap Test To quantify the similarity between the input vector x t and a parameter vector θ, we employ a compact swap test formulation [7, 8]. This test estimates the cosine similarity using quantum state overlaps, encoded as probability amplitudes. Given two 4 real-valued vectors x,θ ∈R n , we define normalized states: φ = ∥x∥ √ Z ,− ∥θ∥ √ Z , Z =∥x∥ 2 +∥θ∥ 2 (3) ψ = 1 √ 2 x 1 ∥x∥ , θ 1 ∥θ∥ ,..., x n ∥x∥ , θ n ∥θ∥ (4) These vectors are then padded (if necessary) to match quantum register lengths as powers of 2. The cosine similarity is estimated using the outcome probability of measuring 0 from the control qubit of a controlled-swap (Fredkin) circuit: cos(θ,x)≈ p 2p 0 − 1, p 0 = Pr(measuring 0)(5) This estimation is executed using repeated quantum circuit simulations over a spec- ified number of shots. The resulting cosine similarity acts as a multifaceted quantum proxy for partial autocorrelation, AR and MA in time series analysis. Algorithm 1 outlines state preparation, and Algorithm 2 details the overall compact swap test computation. Algorithm 1 State Preparation for Compact Swap Test Require: x,θ ∈R n 1: Compute norms: ∥x∥, ∥θ∥ 2: if any norm is zero then 3:return default states [1, 0], [1, 0] 4: end if 5: Compute Z ←∥x∥ 2 +∥θ∥ 2 6: Compute φ← h ∥x∥ √ Z ,− ∥θ∥ √ Z i 7: Initialize ψ ← [] 8: for i = 1 to n do 9:Append to ψ: x i ∥x∥ √ 2 , θ i ∥θ∥ √ 2 10: end for 11: return φ,ψ 5 Algorithm 2 Compact Swap Test Projection Require: Input x,θ ∈R n , number of shots S 1: Compute φ,ψ ← prep-swaptest(x,θ) 2: Normalize φ,ψ to unit ℓ 2 norm 3: Pad ψ to length 2 k for minimal k such that 2 k ≥ len(ψ) 4: Construct quantum registers: one control qubit, one ancilla for φ, and k ancillas for ψ 5: Initialize φ into ancilla qubit, ψ into multi-qubit register 6: Apply Hadamard to control qubit 7: Apply controlled-swap between φ and ψ 8: Apply Hadamard to control qubit 9: Measure control qubit S times to estimate p 0 10: return cos(θ,x)≈ √ 2p 0 − 1 Algorithm 3 Compact Swap Test for Scaled Dot Product Require: Vectors x,θ ∈R n , shots S 1: Compute norms: ∥x∥, ∥θ∥ 2: if ∥x∥ = 0 or ∥θ∥ = 0 then 3:return 0 4: end if 5: Compute state encodings: φ,ψ ← prep-swaptest(x,θ) 6: Normalize: φ← φ/∥φ∥, ψ ← ψ/∥ψ∥ 7: Compute target length: L← 2 ⌈log 2 (len(ψ))⌉ 8: if len(ψ) < L then 9:Pad ψ with zeros to length L 10: end if 11: k ← log 2 (L)▷ Number of qubits for ψ 12: Allocate registers: one control qubit q 1 , 1 +k ancilla qubits q 2 , and 1 classical bit c 13: Initialize φ onto q 2 [0] and ψ onto q 2 [1],...,q 2 [k] 14: Apply: X and H to control qubit q 1 15: Apply: Controlled-swap between q 2 [0] and q 2 [1] 16: Apply: H to q 1 and measure into c 17: Simulate circuit for S shots using QASM backend 18: Obtain measurement probabilities: p 0 ← Pr(0), p 1 ← Pr(1) 19: Compute swap test projection: P swap ← 1− 2p 0 + p 1 20: return ∥x∥·∥θ∥· p max(P swap , 0) 6 |q c ⟩ H • H p 0 |φ⟩ CSWAP |ψ⟩ Fig. 1: Compact swap test estimating similarity between encoded states |φ⟩ and |ψ⟩. The control qubit is prepared in |0⟩, Hadamard → CSWAP → Hadamard, then mea- sured to yield p 0 , from which cosθ swap and the entropy term H(p 0 ) inL AR /L MA are computed. 4.2 Differentiability To determine the optimal differencing order d required to induce stationarity in the time series, we employ a quantum-inspired projection framework in place of tradi- tional unit root tests such as the Augmented Dickey-Fuller (ADF) test [9]. The process begins by constructing a delay matrix from the raw series (Algorithm 4) and gener- ating differenced variants up to a pre-defined maximum order d max (Algorithm 5). For each differencing level d, a pair of parameters (α,γ) are estimated to minimize the mean squared error (MSE) between the differenced target Del d and its predicted value ˆy t , computed via swap-test–based quantum projection. Specifically, we define the projection-based prediction as ˆy t =⟨ψ(x t )|ψ(θ)⟩,(6) where x t = [1,x t ] and θ = [α,γ](7) are encoded into quantum states using amplitude encoding. The objective at each differencing level is to minimize: L d (α,γ) = 1 N X t (Del d (t)− ˆy t ) 2 .(8) The optimization is performed using a classical method such as COBYLA. A con- vergence mechanism is employed: if the learned γ parameter stabilizes over multiple iterations (indicating no further drift), the process halts early. The selected order d ∗ is the smallest d for which the loss is minimized and γ ≈ 0, ensuring that the series is stationary while avoiding unnecessary over-differencing. The complete estimation routine is summarized in Algorithm 6. 7 Algorithm 4 Build Delay Matrix Require: Series y, maximum lag p 1: for i = 1 to p do 2:Create column lag i ← y t−i 3: end for 4: Append target y t return DataFrame of p lag columns and target y t Algorithm 5 Generate Differenced Series Require: Series y, max differencing order d max 1: Set prev← y 2: for i = 1 to d max do 3:Del i ← prev− prev t−1 4:prev← Del i 5: end forreturn DataFrame with Del 1 ,..., Del d max 8 Algorithm 6 Quantum-Inspired Estimation of Differencing Order d Require: Time series y, maximum differencing depth d max , lag order p, convergence threshold ε, patience T , max optimizer iterations I max , primary loss metric L, random seed s 1: Construct delay matrix from y using p lags (Algorithm 4) 2: Generate differenced series up to d max (Algorithm 5) 3: Initialize d← 0, metricslog ←∅ 4: Evaluate d = 0: 5: Let x t ← lag 1 (t), y t ← y(t) 6: Optimize (α 0 ,γ 0 ) to minimize: L 0 (α,γ) = 1 N X t (y t − ˆy t ) 2 7: Compute predictions via swap test projection: ˆy t =⟨ψ([1,x t ])| ψ([α,γ])⟩ 8: Store α 0 , γ 0 , L 0 in metricslog 9: for d = 1 to d max do 10:if Del d or Del d−1 not defined then 11:Break 12:end if 13:Let x t ← Del d−1 (t), y t ← Del d (t) 14:Initialize optimizer with seed s + d 15:Optimize (α d ,γ d ) to minimize: L d (α,γ) = 1 N X t (y t − ˆy t ) 2 16:Compute predictions via swap test projection: ˆy t =⟨ψ([1,x t ])| ψ([α,γ])⟩ 17:Store evaluation metrics in metrics log 18:if last T values of γ d satisfy max− min < ε then 19:Break▷ Early stopping if γ stabilizes 20:end if 21: end for 22: Select d ∗ with lowest L d and γ d ≈ 0 return d ∗ ,α d ∗ ,γ d ∗ The estimation process is governed by several hyperparameters that control con- vergence behavior, optimization stability, and evaluation fidelity. These include the 9 maximum differencing depth d max , maximum number of optimization iterations, con- vergence threshold for γ, a patience parameter that determines early stopping, and a tunable random seed for initialization. Additionally, the primary evaluation metric such as MSE or mean absolute error (MAE) can be selected to align with specific forecasting goals. Table 1: Hyperparameters for Estimating Differencing Order d SymbolNameDescription d max Max differencing depthMaximum number of differencing levels tested. pLag orderNumber of past lags used to construct the delay matrix. εConvergence thresholdThreshold for detecting convergence of γ over recent itera- tions. TPatienceNumber of iterations required to confirm γ stability. I max Max optimizer iterationsMaximum steps allowed per optimization of (α,γ). LLoss metricPrimary evaluation metric (e.g., MSE, MAE) used for model selection. sRandom seedSeed for reproducibility across differencing levels. 4.3 Quantum-Inspired Autocorrelation Function (ACF) The ACF quantifies linear dependence between y t and its lagged values y t−k without conditioning on intermediate lags, and is a primary diagnostic for identifying MA structure in classical ARIMA [1, 10]. We introduce a quantum-inspired estimator that replaces classical inner products by overlap amplitudes computed via a compact swap test [7]. Given a univariate series y t N t=1 (optionally centered and/or differenced), define aligned pairs(y t ,y t−k ) N t=k+1 for lag k. Let ψ(·) denote amplitude encoding of a real scalar into a single-qubit state (or a short register when batching); the swap- test estimate of cosine similarity between encodings is obtained from the control-qubit outcome probability p 0 : [ cosθ t,k = q 2p (t,k) 0 − 1 ,(9) which serves as a quantum proxy for normalized correlation at lag k. Averaging over aligned pairs yields the quantum-inspired ACF value bρ Q k = 1 N − k N X t=k+1 [ cosθ t,k .(10) To enhance robustness, we optionally apply a phase-corrected cosine (as in our PACF module) that blends classical and swap-test angles: cos φ corr t,k = cos θ (t,k) swap + ω θ (t,k) dot − θ (t,k) swap ,(11) 10 and set bρ Q k = 1 N−k P t cos(φ corr t,k ). Significant lags are then selected by thresholding|bρ Q k | using either (i) classical asymptotic bounds±Z/ √ N (with Z = 1.96 for 95% bands), (i) a percentile rule on magnitudes, or (i) a mean+σ cutoff; a fallback band at βτ retains near-significant lags to stabilize downstream MA identification. In classical diagnos- tics, slowly decaying ACF suggests MA order q [1]; here, the quantum-inspired ACF provides the same role but with overlap-based similarity that is naturally compatible with our swap-test losses and VQC refinement. Algorithm 7 Quantum-Inspired ACF Estimation via Compact Swap Test Require: Series y 1:N (optionally centered/differenced), max lag K, shots S, threshold mode T m ∈static, percentile, std, phase weight ω, fallback ratio β Ensure: ACF values bρ Q k K k=1 , selected lags L 1: Initialize list ACF← [ ] 2: for k = 1 to K do 3:if N ≤ k + 1 then 4:Append 0 to ACF; continue 5:end if 6:Form aligned pairs (y t ,y t−k ) N t=k+1 7:Initialize accumulator s← 0, count m← 0 8:for t = k + 1 to N do 9:Encode y t ,y t−k ; run compact swap test with S shots to get [ cosθ t,k 10:Optionally compute classical cosine and apply phase correction to get cos(φ corr t,k ) 11:s← s + cos(φ corr t,k ) or [ cosθ t,k ; m← m + 1 12:end for 13: bρ Q k ← s/m;Append bρ Q k to ACF 14: end for 15: Compute magnitudes a k ← bρ Q k 16: if T m = static then 17: τ ← Z/ √ N 18: else if T m = percentile then 19: τ ← Percentile(a,q) 20: else if T m = std then 21: τ ← mean(a) + σ std(a) 22: else 23:error 24: end if 25: τ f ← β· τ ; L sig ←k : a k ≥ τ; L fb ←k : τ f ≤ a k < τ 26: return ACF =bρ Q k , L←L sig ∪L fb In the classical Box–Jenkins workflow, the ACF guides MA order selection by inspecting damped sine waves and cutoffs [1]. Our quantum-inspired ACF preserves this role but replaces sample correlations with overlap-based similarities that (i) align 11 Table 2: Hyperparameters for Quantum-Inspired ACF Estimation SymbolNameDescription KmaxlagMaximum lag evaluated for ACF. SshotsNumber of swap-test measurements per pair (y t ,y t−k ). T m threshold modeThresholding strategy: static (±Z/ √ N ), percentile, or std. ZthresholdzZ-score for static bounds (typically 1.96 for 95%). qdynamicpercentilePercentile for magnitude cutoff (e.g., 90th). σstdmultiplierMultiplier for std-based bound (mean + σ·std). ωphaseweightWeight for phase-corrected cosine blend (optional). βfallbackratioFallback fraction of the primary threshold (0 < β < 1). natively with our swap-test losses used later in MA training, and (i) integrate a principled uncertainty control through shot noise and optional entropy regularization in downstream objectives. Together with the quantum-inspired PACF, this yields a consistent diagnostic pair for proposing candidate (p,q) before VQC-based refinement. 4.4 Quantum-Inspired PACF Estimation. The PACF is a fundamental diagnostic tool in time series analysis used to determine the order p of the autoregressive (AR) model by measuring the correlation between y t and y t−k , after eliminating the influence of intermediate lags [1]. In this work, we extend PACF estimation through a quantum-inspired approach that leverages projection-based inner products via the compact swap test, a method to quantify the similarity between the input vector x t and a parameter vector θ in a quantum-inspired PACF estimation [7, 8] Algorithms 1, 2. For a given lag k, the projection amplitude is computed as ˆρ k = 1 N − k N X t=k+1 ⟨ψ(y t )|ψ(y t−k )⟩,(12) where ψ(·) denotes a quantum state encoding of the input scalar into qubit ampli- tude representations [11]. The swap test is repeated over N−k aligned pairs to generate a mean partial correlation value at lag k. After computing PACF values across all lags up to a maximum K, a thresholding mechanism is applied to select significant lags. We implement three strategies: • a static threshold τ = 1.96 √ N based on classical asymptotic confidence intervals [10] • a dynamic threshold based on the q-th percentile of PACF magnitudes • a standard deviation threshold defined as μ + σ, where μ and σ are the mean and standard deviation of the PACF magnitudes To enhance robustness, a fallback mechanism selects lags with projections between a fallback threshold τ f = β · τ and the primary threshold τ , ensuring the retention of moderately informative lags. This framework preserves the interpretability of classical 12 PACF while embedding quantum-aligned projection logic into lag selection, a tech- nique not commonly adopted in the time series literature, thereby introducing a novel avenue for hybrid quantum-classical model design. Algorithm 8 Quantum-Inspired PACF Estimation Using Swap Test Require: Time series y, maximum lag K, swap test function SwapTest, number of shots S, threshold mode T m , fallback flag f 1: Initialize PACF← [] 2: for k = 1 to K do 3:if length(y)≤ k + 1 then 4:Append 0 to PACF and continue 5:end if 6:Extract y t ←y k+1 ,...,y N , y t−k ←y 1 ,...,y N−k 7:for each t do 8:Compute projection: ρ (k) t =⟨ψ(y t )|ψ(y t−k )⟩≈ SwapTest(y t ,y t−k ,S) 9:end for 10:Compute ˆρ k = 1 N−k P t ρ (k) t 11:Append ˆρ k to PACF 12: end for 13: Compute threshold τ : 14: if T m = “static” then 15: τ ← Z √ N 16: else if T m = “percentile” then 17: τ ← top-q percentile of |PACF| 18: else if T m = “std” then 19: τ ← μ + σ 20: end if 21: Define fallback threshold τ f ← β· τ 22: Identify significant lags: L s =k || ˆρ k |≥ τ, L f =k | τ f ≤| ˆρ k | < τ 23: Combine: L final ←L s ∪L f if f = True return PACF, L final , τ 4.5 AutoRegressive Order p In ARIMA modeling, the autoregressive (AR) order p governs how many past values of a time series y t are used to predict its current value. Classical AR models rely on 13 Table 3: Hyperparameters for Quantum-Inspired PACF Estimation SymbolNameDescription KmaxlagMaximum lag to compute PACF for. SshotsNumber of measurement shots used in the swap test circuit. T m thresholdmodeStrategy to determine the threshold: “static”, “percentile”, or “std”. ZthresholdnumeratorZ-score for static thresholding (typically 1.96 for 95%). qdynamicpercentilePercentile threshold (e.g., 90th percentile) for projection magni- tude. σstdmultiplierNumber of standard deviations above mean (used in “std” mode). βfallbackratioProportion of τ used as the fallback threshold. fenablefallbackBoolean flag to include fallback lags in final selection. linear regression, using dot-product formulations over lagged vectors. In our quantum- inspired extension, we replace the classical regression-based fitting with a loss function that combines classical and quantum elements: phase-corrected cosine similarity from compact swap tests, entropy-based uncertainty, and misalignment penalties. The pro- cess is formulated as follows. For a lag vector x t and coefficient vector b at time t, the classical prediction is: ˆy (dot) t = x ⊤ t b(13) We compute the classical cosine similarity as: cosθ dot = x ⊤ t ˆ b ∥x t ∥ ,where ˆ b = b ∥b∥ (14) The quantum-inspired variant estimates the cosine angle using a simulated swap test: cosθ swap = SWAP(x t , b)(15) This estimate is phase-corrected based on the discrepancy between θ dot = arccos(cosθ dot ) and θ swap = arccos(cosθ swap ), resulting in the corrected prediction: ˆy (quantum) t =∥x t ∥·∥b∥· cos (θ swap + λ phase (θ dot − θ swap ))(16) The total AR loss integrates prediction error, cosine similarity misalignment, and binary entropy uncertainty from the swap test probability: L AR = X t y t − ˆy (quantum) t 2 + λ cos X t (cosθ dot − cosθ swap ) 2 + λ ent X t H(cosθ swap ) (17) 14 where the binary entropy is computed as: H(p) =−p log 2 p− (1− p) log 2 (1− p)(18) The following algorithms compute the total loss L AR for a given p and coefficient vector b: Algorithm 9 Phase-Corrected Cosine Calculation Require: Classical cosine cosθ dot , Quantum cosine cosθ swap , Phase correction weight ω Ensure: Corrected cosine value and intermediate angles 1: Clip classical cosine: ̃c dot ← clip(cosθ dot ,−1.0, 1.0) 2: Clip quantum cosine: ̃c swap ← clip(cosθ swap ,−1.0, 1.0) 3: Compute angles: θ dot ← arccos( ̃c dot ), θ swap ← arccos( ̃c swap ) 4: Phase difference: ∆θ ← θ dot − θ swap 5: Corrected angle: θ corr ← θ swap + ω· ∆θ 6: Corrected cosine: cosθ corr ← cos(θ corr ) 7: return cosθ corr ,θ dot ,θ swap , ∆θ,θ corr Algorithm 10 Quantum-Inspired AR Loss Estimation Require: Lagged features x t , targets y t , candidate AR coefficients b, swap test function, hyperparameters λ cos ,λ ent ,λ phase Ensure: Total quantum-inspired AR loss L AR 1: Normalize b to obtain unit vector b unit 2: Pad b unit and x t to power-of-2 length 3: for each time step t do 4:Compute classical dot product: ˆy (dot) t ← x ⊤ t b 5:Compute classical cosine: cosθ dot ← x ⊤ t b unit ∥x t ∥ 6:Estimate cosθ swap via compact swap test 7:Compute angular correction: ∆θ ← θ dot − θ swap 8:Corrected angle: θ corr ← θ swap + λ phase · ∆θ 9:Quantum prediction: ˆy (quantum) t ←∥x t ∥·∥b∥· cos(θ corr ) 10:Entropy: H t ←−p 0 log 2 p 0 − (1− p 0 ) log 2 (1− p 0 ) where p 0 = 1− cos 2 θ swap 11:Accumulate: L += (y t − ˆy (quantum) t ) 2 + λ cos (cosθ dot − cosθ swap ) 2 + λ ent H t 12: end for 13: return L 15 Table 4: Hyperparameters for Quantum-Inspired AR Loss Estimation SymbolNameDescription SshotsNumber of swap test measurement shots used for cosine similarity estimation. λ cos lambdacosinealignPenalty weight applied to alignment difference between classical dot and swap-based cosine. λ ent lambdaentropyPenalty weight for quantum entropy derived from swap test prob- ability. ωphaseweightScaling factor to correct the phase discrepancy between classical and quantum cosine. σnoisestdStandard deviation of added Gaussian noise to simulate quantum fluctuations. PpaddeddimSize of padded vector, adjusted to nearest power-of-two for swap test register. fsavecsvBoolean flag to save log metrics and loss values for each AR trial. 4.5.1 AR Order Estimation via VQC Refinement In classical time series analysis, the autoregressive order p is typically selected by scanning information criteria such as AIC or BIC, or by grid search guided by PACF [1]. In our quantum-inspired framework, we preserve the PACF-guided candidate set P, but replace global search methods with a VQC refinement of the AR coefficients for each candidate p. For each p∈P, we proceed as follows: 1. Difference the series using the identified order d ⋆ ; 2. Build the AR design matrix X (p) from lagged values; 3. Initialize the coefficient vector b (0) via ordinary least squares (OLS); 4. Refine b by minimizing the quantum-inspired AR loss L AR (cf. Eq. (17)) using a depth-controlled VQC ansatz and a classical derivative-free optimizer such as COBYLA [4, 12]. The VQC offers an expressive, hardware-compatible parameterization while retain- ing PACF interpretability. OLS initialization provides a warm start that reduces susceptibility to barren plateaus and poor local minima [13]. Let y ∈R T denote the differenced series. For each p∈P, we form X (p) ∈R (T−p)×p and solve: b (0) = arg min b ∥y− X (p) b∥ 2 2 .(19) The initialized coefficients are embedded into the VQC as R y (b j ) rotations. The VQC-based AR loss is defined as: L VQC AR (β) = X t y t − ˆy (quantum) t (β) 2 + λ cos X t (cosθ dot,t − cosθ swap,t ) 2 + λ ent X t H(p 0,t ), (20) where the phase-corrected prediction is: ˆy (quantum) t (β) = ∥x t ∥b(β)∥ cos(θ swap,t + ω (θ dot,t − θ swap,t )),(21) 16 with: cosθ dot,t = x ⊤ t b b(β) ∥x t ∥ ,(22) b b(β) = b(β) ∥b(β)∥ ,(23) p 0,t = 1− cos 2 θ swap,t ,(24) H(p) = −p log 2 p− (1− p) log 2 (1− p).(25) The optimal AR order is chosen as: p ⋆ = arg min p∈P L VQC AR b ⋆ (p) ,(26) where b ⋆ (p) are the coefficients obtained after VQC refinement. Algorithm 11 VQC-Based AR Order Estimation Require: Differenced series y, best differencing order d ⋆ , candidate setP, number of VQC layers r, maximum iterations T max Ensure: Best AR order p ⋆ , optimized coefficients b ⋆ 1: for each p∈P do 2:Build delay matrix X (p) from y 3:Initialize b (0) via OLS 4:Embed b (0) into VQC rotation angles 5:Minimize L VQC AR (β) via COBYLA for T max iterations 6:Store (p, b ⋆ (p) ,L VQC AR (b ⋆ (p) )) 7: end for 8: p ⋆ ← arg min p∈P L VQC AR (b ⋆ (p) ) 9: return p ⋆ , b ⋆ (p ⋆ ) Table 5: Hyperparameters for VQC-Based AR Order Estimation SymbolNameDescription PcandidatelagsPACF-identified set of candidate AR orders. d ⋆ bestdOptimal differencing order from quantum ADF tests. b (0) b initOLS-initialized AR coefficients. rvqclayersNumber of VQC entangling layers. SshotsNumber of measurement shots in swap tests. T max maxiterMaximum iterations for VQC optimization. λ cos ,λ ent ,ωloss weightsRegularization and phase correction weights for the loss function. 17 q 1 R y b (0) 1 R y (β (1) 1 ) • R y (β (2) 1 ) • ˆy t =⟨O⟩ q 2 R y b (0) 2 R y (β (1) 2 )R y (β (2) 2 ) . . . . . . . . . . . . q p R y b (0) p R y (β (1) p )R y (β (2) p ) Fig. 2: VQC for AR order estimation with p qubits (one per lag). The first rotation layer uses OLS initialization R y (b (0) i ); then r=2 variational layers (shown) alternate trainable R y (β (ℓ) i ) with a linear entangling pattern (CNOT ladder). Readout yields ˆy t = ⟨O⟩, and parameters are trained to minimize L AR (prediction error + cosine misalignment + entropy). 4.5.2 Weak Lag Refinement for Extended AR Order (VQC-Based) Even after selecting an optimal AR order p ⋆ via PACF-guided, VQC-refined esti- mation, higher-order terms can carry weak but non-negligible predictive signal under sparsity. We therefore extend the anchor model by adding a small set of weak lags and refining their coefficients with a variational procedure while keeping the anchor coefficients fixed. Selecting weak-lag initializations: Let b ⋆ be the AR coefficients for order p ⋆ , and let S denote the set of VQC-refined coefficient vectors collected across tested orders in the summary table. We extract a pool of candidate magnitudes from S, remove those already present among the anchor magnitudes |b ⋆ |, and keep the top k smallest remaining values as weak-lag initializations: w = w 1 ,...,w k , p ′ = p ⋆ + k. This yields an extended AR(p ′ ) design with anchor part (fixed) and k weak-lag coefficients to be refined. Objective with adaptive penalties: Let b weak ∈R k denote the weak-lag coefficients and b full = [ b ⋆ ; b weak ] the concate- nated vector. Using the quantum-inspired AR loss L AR (·) from Eq. (17), we optimize only b weak by minimizing L weak (b weak ) = L AR b full + λ dev b weak − w 2 2 + λ mag b weak 1 ,(27) where the deviation term encourages proximity to the weak-lag initializations w and the magnitude term promotes sparsity. In practice, we scale the penalty weights 18 relative to a baseline loss L base (the best loss in the summary table) via λ dev = 10 −3 L base , λ mag = 5× 10 −4 L base .(28) VQC parameterization: The weak-lag block is encoded by a depth-r VQC ansatz (e.g., R y layers with entangle- ment). The rotation angles correspond to entries of b weak , while the anchor block b ⋆ remains fixed. A classical derivative-free solver (COBYLA) refines b weak by minimizing L weak , consistently with our hybrid training loop. Algorithm 12 VQC-Based Weak Lag Extension and Refinement Require: Time series y, differencing order d ⋆ , best AR order p ⋆ , anchor coefficients b ⋆ , loss summary Summary, weak count k, VQC depth r, max iterations T max Ensure: Extended order p ′ , refined coefficients b ⋆ (p ′ ) , refined loss L weak 1: Extract candidate magnitudes from Summary; remove magnitudes in |b ⋆ |; keep k smallest as w =w 1 ,...,w k 2: Set p ′ ← p ⋆ + k 3: Difference y by d ⋆ and build AR design X (p ′ ) and targets 4: Initialize b (0) weak ← w and form b (0) full ← [ b ⋆ ; b (0) weak ] 5: Build a depth-r VQC with k trainable angles for b weak 6: Define L weak as in Eq. (27); fix b ⋆ , optimize over b weak 7: Minimize L weak using COBYLA (budget T max ), obtain b ⋆ weak 8: Form b ⋆ (p ′ ) ← [ b ⋆ ; b ⋆ weak ] and evaluate L weak 9: return p ′ , b ⋆ (p ′ ) , L weak Table 6: Hyperparameters for VQC-Based Weak Lag Refinement SymbolNameDescription p ⋆ bestpAR order chosen by PACF + VQC refinement. b ⋆ anchorcoeffsFixed anchor coefficients for AR(p ⋆ ). kweakcountNumber of weak lags added to form p ′ = p ⋆ + k. wweakinitsWeak-lag initialization magnitudes (from summary, excluding |b ⋆ |). rvqcdepthNumber of entangling layers in the weak-lag VQC block. T max maxiterMax COBYLA iterations for minimizing L weak . λ dev deviation penaltyWeight on ∥b weak − w∥ 2 2 (scaled by L base ). λ mag magnitudepenaltyWeight on ∥b weak ∥ 1 (scaled by L base ). 19 Anchor block (frozen) Weak-lag block (trainable) R y (β (1) w,1 ) • R y (β (2) w,1 ) • R y (β (1) w,2 ) R y (β (2) w,2 ) R y (β (1) w,3 )R y (β (2) w,3 ) Fig. 3: VQC for weak-lag refinement with p ⋆ frozen anchors (top block) and k trainable weak-lag qubits (bottom block). Anchors are initialized as R y b (p ⋆ ) i and kept fixed (shown as a grouped “frozen” block). Only weak-lag qubits carry trainable R y (β (ℓ) w,j ) and are entangled (CNOT ladder restricted to the weak block). Training minimizes L weak =L AR + λ dev ∥b weak − w∥ 2 2 + λ mag ∥b weak ∥ 1 . 4.5.3 Progressive Weak Lag Refinement via VQC After identifying a viable extended lag order p ′ , we further decompose its structure through a progressive refinement strategy. Starting from the best anchor coefficients b (p ⋆ ) obtained from PACF-guided VQC refinement, we incrementally add one weak lag at a time and re-optimize the extended coefficient vector. This stepwise inclusion facilitates interpretability of higher-order autoregressive models while minimizing the risk of overfitting by allowing early termination if no improvement is observed. Let w = w 1 ,w 2 ,...,w k denote the initialization values for the k weak lags, ranked by their absolute magnitude. For each p = p ⋆ + 1 to p ⋆ + k, we extend the model to include the first j = p− p ⋆ weak lags, forming an initial coefficient vector b (p) init = h b (p ⋆ ) ,w 1 ,...,w j i . We employ a VQC-based refinement over the weak-lag block, keeping the anchor coefficients fixed. The weak-lag optimization problem is formulated as: L (p) weak = L AR b (p) + λ dev b (p) weak − w 1:j 2 2 + λ mag b (p) weak 1 ,(29) where: 20 • L AR is the quantum-inspired AR loss defined in Eq. (17), • b (p) weak denotes the coefficients of the weak-lag block, • λ dev penalizes deviation from the weak-lag initializations, • λ mag encourages sparsity. The penalties are adaptively scaled by the baseline loss L base =L AR (b (p ⋆ ) ): λ dev = 10 −3 L base , λ mag = 5× 10 −4 L base .(30) For each refinement step, the weak-lag parameters are embedded into a j-qubit VQC with R y rotations and r entangling layers, and optimized using COBYLA. The process continues progressively until the loss improvement ∆ (p) =L (p) weak −L base falls below a predefined threshold, at which point refinement halts. Algorithm 13 Progressive Weak Lag Refinement via VQC Require: Differenced series y; best differencing order d ⋆ ; best AR order p ⋆ ; anchor coefficients b (p ⋆ ) ; weak-lag inits w; VQC depth r; max iterations T max ; threshold τ Ensure: Sequence (p, b (p) ,L (p) ) 1: L base ←L AR (b (p ⋆ ) ) 2: for j = 1 to |w| do 3: p← p ⋆ + j 4:Initialize b (p) init ← [b (p ⋆ ) ,w 1 ,...,w j ] 5:Build a j-qubit VQC (depth r) for the weak-lag block 6:Minimize L (p) (Eq. (29)) via COBYLA for T max iterations 7:Record (p, b (p) ,L (p) weak ) 8:if L (p) weak −L base > τ then 9:break 10:end if 11: end for Table 7: Hyperparameters for Progressive VQC-Based Weak Lag Refinement SymbolNameDescription wweaklagsOrdered list of weak-lag initialization values. rvqcdepthNumber of entangling layers in weak-lag VQC. T max maxiterMaximum COBYLA iterations per refinement step. λ dev ,λ mag penaltyweightsAdaptive regularization weights based on L base . 21 4.6 Residual Evaluation for Quantum-Inspired AR Models Once the optimal autoregressive order p and the associated coefficient vector b are determined via quantum-inspired loss minimization, the next step is to compute resid- uals. These residuals represent the error between predicted and observed values of the differenced time series ∆ d y t and are essential for both evaluating AR model accuracy and informing subsequent MA model design. For each lag order p, a delay matrix X t is constructed from p prior values of the differenced series ∆ d y t . The quantum-inspired prediction for each time step is formulated as: ˆy t = ∥b∥·∥X t ∥· cos (θ t )(31) cos(θ t ) = X t · ˆ b ∥X t ∥ (32) ˆ b = b ∥b∥ (33) The residuals are then computed as: ε t = ∆ d y t − ˆy t (34) From the residual sequence ε t , we calculate the mean μ ε , standard deviation σ ε , and preserve the full residual vector for each p in a structured summary. This residual analysis serves as a bridge to MA modeling, where these residuals are treated as noise inputs for estimating MA coefficients. The stepwise procedure is detailed in Algorithm 14. 22 Algorithm 14 Residual Evaluation for Quantum-Inspired AR Models Require: Differenced time series ∆ d y, optimal differencing order d, set of AR orders p i and coefficients b i Ensure: Residual statistics for each AR model 1: Initialize residual log R← [ ] 2: for each row i in AR results do 3: p← p i , b← b i 4:Construct delay matrix X using p lags from ∆ d y 5:Set target vector y target ← ∆ d y p: 6:Normalize b to get unit vector ˆ b 7:for each row x t in X do 8:Compute norm ∥x t ∥ 9:Compute cosine projection: cos(θ t )← x t · ˆ b ∥x t ∥ 10:Compute prediction: ˆy t =∥b∥·∥x t ∥· cos(θ t ) 11:Compute residual: ε t = y t − ˆy t 12:end for 13:Compute μ ε ,σ ε ← mean and std of residuals 14:Append (p,μ ε ,σ ε ,ε t ) to R 15: end for 16: return Residual summary table R 4.7 Quantum-Inspired Moving Average Modeling In classical time series modeling, MA components model the serial correlation in noise or residual errors, capturing short-term dependencies not explained by autoregressive terms. In this work, we introduce a quantum-inspired MA modeling framework that leverages quantum similarity estimation via swap tests, phase-corrected projections, and entropy-aware regularization. Given a time series y t T t=1 and AR residuals ε t derived from the autoregressive model, our approach performs a hybrid modeling of residuals as an MA(q) process using variational optimization informed by quantum projection. Specifically, we define a trainable MA parameter vector θ ∈R q and seek to minimize a quantum-regularized loss function over candidate values of q. The overall MA modeling framework includes the following steps: • Loss computation: Uses quantum projection and phase correction to estimate the predictive error for a given θ. • MA order estimation: Explores various values of q and selects the best one via minimal total loss. • Model training: Applies constrained optimization to train the MA model for the selected q. • ARMA extension: Iterates across AR models and fits MA components to generate the final ARMA model summary. Each of these components is described in detail below. 23 4.7.1 Quantum MA Loss Estimation To estimate the quality of MA coefficients θ, we define a quantum-regularized loss function L MA that integrates classical residual prediction with quantum projections derived from swap tests [14, 15]. This formulation enables us to model the noise corre- lation in residuals of an AR process while incorporating quantum-inspired similarity estimates. We define the delayed residual vector at time t as: ε t = [ε t−1 ,ε t−2 ,...,ε t−q ](35) and the predicted target as: ˆy t =∥θ∥·∥ε t ∥· cos(φ corrected t )(36) Here, φ corrected t is a phase-adjusted angle between the normalized θ and ε t , combining classical and quantum projections: cos(φ corrected t ) = cos (θ swap + ω(θ dot − θ swap ))(37) where: • θ dot = arccos(⟨θ unit ,ε t,unit ⟩) is the classical angle, • θ swap is estimated via the compact swap test, • ω is a tunable phase correction weight. The total quantum MA loss function is given by: L MA = X t (y t − ˆy t ) 2 + λ cos X t cos(φ t )− cos(φ corrected t ) 2 + λ ent X t H t + λ L2 ∥θ∥ 2 (38) where H t is the entropy of the swap-test-derived distribution: H t =− X i p i log 2 (p i ),with p 0 = cos 2 (φ swap t ), p 1 = 1− p 0 (39) 24 Algorithm 15 Quantum-Inspired MA Loss Estimation Require: MA coefficients θ, residual matrix E, target series y t , trial ID, reg- ularization weights λ cos ,λ ent ,λ L2 , phase weight ω, number of quantum shots S Ensure: Quantum-regularized loss L MA 1: Normalize θ to unit vector θ unit 2: Pad all vectors to nearest power-of-2 length for swap test 3: for each time step t do 4:Retrieve residual window ε t 5:Compute classical cosine: cos(θ dot ) via dot product 6:Estimate quantum cosine: cos(θ swap ) via compact swap test 7:Apply phase correction using Eq. (3) to get cos(φ corrected t ) 8:Predict: ˆy t =∥θ∥·∥ε t ∥· cos(φ corrected t ) 9:Compute squared error, cosine deviation penalty, and entropy H t 10:Accumulate total loss components 11: end for 12: Add λ L2 ·∥θ∥ 2 regularization 13: return L MA Table 8: Hyperparameters for Quantum MA Loss SymbolNameDescription λ cos lambdacosinealignPenalty for misalignment between classical dot-product and quan- tum cosine. λ ent lambdaentropyWeight for entropy regularization based on swap test measure- ments. ωphase weightPhase correction factor for angular adjustment. SshotsNumber of measurements used in the swap test simulation. λ L2 lambdal2ℓ 2 regularization weight on the magnitude of θ. 4.7.2 MA Order Estimation via Quantum Residual Modeling To determine the optimal MA order q, we evaluate the quantum-inspired MA loss L MA across a candidate set Q = 1, 2,...,q max . For each value of q, the model fits MA coefficients θ (q) to the delayed residuals using a variational optimizer guided by quantum projection loss. The best order q ⋆ is selected based on the minimal loss: q ⋆ = arg min q∈Q L (q) MA (40) 25 Each candidate q generates a delay matrix from residuals, against which the pre- diction quality is assessed using the previously defined loss function. Initialization is done by sampling the MA coefficients from a standard Gaussian distribution: θ (0) ∼N (0,I q )(41) The loss function L MA incorporates projection-based cosine similarity using the swap test [7], phase correction to improve alignment between classical and quantum cosine [14], and entropy-based regularization to manage uncertainty from quantum measurements [15]. Optimization of each θ (q) is performed using classical constrained routines such as COBYLA [12], consistent with variational quantum algorithms in hybrid settings [4]. Algorithm 16 Quantum MA Order Estimation Require: Residuals ε t , target series y t , candidate range Q = [q min ,q max ], loss function L MA Ensure: Optimal MA order q ⋆ 1: Initialize empty summary list 2: for each q ∈Q do 3:Construct delay matrix E q using q lags 4:Initialize MA vector θ (q) ∼N (0,I q ) 5:Compute loss: L (q) MA ← QuantumMAProjectionLoss(θ (q) , E q ,y t ) 6:Store (q,L (q) MA ) in summary 7: end for 8: Select q ⋆ ← arg minL (q) MA 9: return q ⋆ Table 9: Hyperparameters for MA Order Estimation SymbolNameDescription q min ,q max qrangeRange of candidate MA orders tested during estimation. θ init thetainitRandom initialization of coefficient vector for each trial q. L MA lossfunctionQuantum-inspired MA loss computed using swap test projections. 4.7.3 MA Coefficient Optimization via VQC Refinement After selecting the optimal MA order q ⋆ via the QACF-gated procedure, we perform a final parameter estimation by minimizing the quantum-augmented loss L MA (θ) with 26 respect to the coefficient vector θ ∈R q ⋆ . The loss retains the original structure in Equation (38), comprising: 1. Prediction error computed using phase-corrected cosine similarity, 2. Cosine alignment penalty between dot-product and swap-test derived similarities, 3. Entropy penalty on swap-test projections, 4. ℓ 2 regularization for coefficient stability. The optimization is formulated as a norm-constrained problem: min θ L MA (θ) s.t. ∥θ∥ 2 ≤ τ,(42) where τ is an upper bound to mitigate overfitting and improve generalization. VQC-based parameterization: Instead of directly optimizing θ in Euclidean space, we embed its components as rotation parameters in a shallow VQC denoted buildvqcma(q ⋆ ,r), where r is the number of entangling repetitions. This provides a hardware-compatible representation while preserving a one-to-one correspondence between the parameter vector and the circuit’s R y rotation angles. The optimiza- tion is then carried out in a hybrid loop: the VQC generates quantum projections for loss evaluation, and a classical derivative-free solver (COBYLA [12]) updates the parameters to minimize Eq. (42). Initialization. We warm-start θ using conditional least squares (Hannan– Rissanen). First, from the residual stream ε t we form the delayed error matrix E (q ⋆ ) = [ε t−1 · ε t−q ⋆ ] aligned to the MA target vector y MA . We then solve θ (0) = arg min θ y MA − E (q ⋆ ) θ 2 2 ,(43) and clip the entries elementwise to [−1, 1] for numerical stability. If CLS initialization fails (e.g., ill-conditioned design or insufficient data), we fall back to θ (0) ∼U (−1, 1). Post-training evaluation: Upon convergence, we compute the empirical stan- dard deviation of the residuals to quantify the remaining volatility in the MA model: σ MA = v u u t 1 T T X t=1 (y t − ˆy t ) 2 ,(44) where ˆy t is the VQC-predicted output using the optimized θ ⋆ . This measure, together with the convergence loss, serves as a post-training quality indicator. 27 Algorithm 17 VQC-Refined MA Training for q ⋆ Require: Residuals ε t , AR target y t , optimal order q ⋆ , VQC depth r Ensure: Optimized coefficients θ ⋆ , MA noise standard deviation σ MA 1: Construct delayed error matrix E (q ⋆ ) from ε t 2: Initialize θ (0) via OLS; fallback to U (−1, 1) 3: Embed θ into buildvqcma(q ⋆ ,r) 4: Minimize L MA (θ) via COBYLA with ∥θ∥ 2 ≤ τ 5: Compute ε t = y t − ˆy t and σ MA = q 1 T P t ε 2 t 6: return θ ⋆ , σ MA Table 10: Hyperparameters for VQC-Refined MA Training SymbolNameDescription θ (0) thetainitInitial coefficient vector from OLS (clipped) or U (−1, 1). τnormconstraintNorm bound on θ to control overfitting. rvqcrepsNumber of entangling layers in the VQC ansatz. T max maxiterMaximum iterations for COBYLA. SshotsNumber of swap-test shots for loss evaluation. λ cos ,λ ent ,λ L2 ,ωlossweightsLoss weights and phase-correction factor (as in Eq. 38). q 0 R y (θ 0 ) • · q 1 R y (θ 1 ) • · q 2 R y (θ 2 ) · . . . q q ⋆ −2 R y (θ q ⋆ −2 ) • · q q ⋆ −1 R y (θ q ⋆ −1 ) · Fig. 4: Schematic VQC ansatz for MA coefficient refinement buildvqcma(q ⋆ ,r). Each qubit q i receives a single-qubit R y (θ i ) rotation, followed by an entangling CNOT chain (pattern shown) and a block repeated r times (indicated by · ). The resulting quantum state is used to compute L MA . 28 4.7.4 ARMA Finalization Across All AR Models Finally, the full ARMA model is built by applying the MA estimation procedure to each AR(p) model already selected through quantum optimization. The final DataFrame summarizes (p,d,q) combinations and their associated AR and MA parameters, allowing for direct comparison of ARMA configurations. Algorithm 18 Run MA Estimation for All AR Models Require: Differenced data y, AR model summary (p,b), residual matrix Ensure: Final DataFrame of ARMA (p,d,q) models 1: for each AR model with order p do 2:Retrieve AR residuals 3:Apply MA estimation and training 4:Store (p,d,q,b,θ) 5: end for 6: return DataFrame of ARMA models Table 11: Output Structure for Final ARMA Summary ColumnTypeDescription pIntegerAR order dIntegerDifferencing order qIntegerMA order (quantum-optimized) bVectorAR coefficients θVectorMA coefficients σ AR FloatNoise std from AR residuals σ MA FloatEstimated MA noise std 5 Results: Evaluation Protocol Across Datasets We evaluated Quantum-ARIMA (QARIMA) against a classical pmdarima baseline on five public time series: Sunspots, Mauna Loa CO 2 , Australian Woollen Yarn Produc- tion, Australian Beer Production, and Sydney 2024 Weather. For each dataset, we followed the Methods pipeline: (i) identical preprocessing and differencing, (i) gener- ation of candidate QARIMA(p,d,q) models via the manuscript procedure, and (i) training of a matched classical ARIMA comparator on the same training data. Model quality was assessed out-of-sample (OOS) on a held-out segment defined per dataset. We report error metrics (MSE, MAPE) for every candidate model and then conduct Diebold-Mariano (DM) tests comparing each QARIMA against the classical baseline under two loss functions (MSE and MAE). 29 5.1 Error-Based Evaluation For every dataset we report out-of-sample (OOS) accuracy using two standard ARIMA metrics: Mean Squared Error (MSE) and Mean Absolute Percentage Error (MAPE). MSE reflects point-forecast fidelity and penalizes large deviations, while MAPE pro- vides a scale-free view that is comparable across series with different magnitudes. For each QARIMA candidate we list its OOS value alongside the matched classical baseline and, when useful, the absolute improvement ∆MSE = MSE classical − MSE quantum ,∆MAPE = MAPE classical − MAPE quantum . Individual dataset subsections interpret these numbers in context (trend-like vs. sea- sonal series, short vs. long OOS windows) and highlight cases where both metrics move in the same favourable direction. 5.2 Statistical Comparison To determine whether the observed error gaps are more than numerical fluctuations, we run Diebold-Mariano (DM) tests [16] for each quantum–classical pair on the same OOS segment. We evaluate DM under two loss differentials: squared-error (MSE-style) and absolute-error (MAE-style). For every comparison we report the DM statistic, its p- value, and mark results that satisfy p≤ α (typically α = 0.05). Alongside significance, we also show the corresponding classical and quantum mean losses and their difference ∆ = loss classical − loss quantum , which acts as an effect-size indicator. The per-dataset result sections that follow use this common procedure to comment on when a QARIMA variant is not only better in value (MSE/MAPE) but also statistically supported by DM. Interpretation rule. We consider a QARIMA model to reliably outperform the classical baseline on a dataset when it (i) improves OOS MSE/MAPE and (i) achieves DM significance (p ≤ 0.05) under at least one loss (preferably both). When multiple QARIMA speci- fications are significant, we prioritize those with the largest positive mean-loss deltas and consistent gains across both error metrics. In subsequent subsections, we will present the results achieved by each dataset and our interpretation of QARIMA as per the results. 5.3 Sunspots dataset Weusetheclassicannualmeansunspotcountsprovidedin statsmodels.datasets.sunspots (column SUNACTIVITY), spanning 1700-2008. These counts proxy solar magnetic activity and exhibit the well-known ∼11-year Schwabe cycle, but with quasi-periodic behavior: cycle length and amplitude drift over time, peaks are asymmetric, and multi-decadal envelopes (amplitude/phase mod- ulation) are common. Statistically, the series is nonstationary in level and nonlinear 30 Table 12: Sunspots Classical Vs Quantum MSE MAPE Sunspots Classical Vs Quantum OOS MSE MAPE ModelNMSEMAPE Classical (pmdarima) (2, 0, 0) seasonal=None1282181.5891.661 Quantum (pdq=(3, 1, 1))1282146.9261.790 Quantum (pdq=(10, 1, 1))1282160.9981.534 Quantum (pdq=(7, 1, 3))1282174.6221.629 Quantum (pdq=(4, 1, 1))1282179.4051.705 Quantum (pdq=(9, 1, 3))1282182.9781.601 Quantum (pdq=(8, 1, 3))1282184.0061.580 Quantum (pdq=(6, 1, 3))1282210.4341.349 Quantum (pdq=(5, 1, 3))1282256.7881.254 in its dynamics; simple linear AR models often underfit long-range dependence and changing cycle shape, motivating differencing and richer lag structure. Sunspots have total of 308 annual data points out of which we have use 181 data points for training and 128 data points for Out of Sample test. The generated QARIMA models are evaluated against the classical pmdarima(2,0,0). The OOS results of MSE , MAPE and DM tests are presented in Table 12. Sunspots OOS performance assessment. MSE / MAPE comparison. All QARIMA models were benchmarked against the classical non-seasonal pmdarima ARIMA(2, 0, 0) on the 128-step OOS window. The results are shown in Table 12 and in the bar charts in Fig. 5. Several quantum specifi- cations clearly match or beat the classical baseline on MSE, most notably Q(10,1,1), Q(9,1,3), Q(8,1,3), Q(7,1,3), Q(6,1,3), and Q(5,1,3), which all reduce OOS MSE rel- ative to ARIMA(2, 0, 0) while keeping MAPE in a comparable range. Two models, Q(6,1,3) and Q(5,1,3), show the largest MAPE improvement (down to ≈ 1.35 and ≈ 1.25), indicating that the quantum variants can better track relative deviations in the tail of the 128-step sequence. By contrast, Q(3,1,1) and, to a lesser extent, Q(4,1,1), stay close to or slightly above the classical MSE/MAPE and do not contribute to the best-performing group. DM analysis. To verify that the observed gains are not due to sampling variabil- ity, we ran Diebold–Mariano tests between each QARIMA forecast and the classical ARIMA(2, 0, 0) forecast. On the MSE loss, most quantum models (Q(10,1,1), Q(9,1,3), Q(8,1,3), Q(7,1,3), Q(6,1,3), Q(5,1,3), and Q(4,1,1)) fall below the α = 0.05 line in the DM panels (Fig. 6), confirming that their lower mean loss is statistically significant over the 128 OOS points. The only consistent non-significant case is Q(3,1,1), whose error sequence is too close to the classical one. On the MAE/DM side (Fig. 7) the pattern is similar but slightly weaker: 6 of the 8 quantum models clear the α = 0.05 threshold, led again by the higher-order quantum specifications; Q(3,1,1) and occa- sionally Q(4,1,1) remain above the line. Full DM statistics (DM stat, p-value, classical and quantum mean losses, and ∆) are provided in the appendix for reproducibility. 31 Fig. 5: Sunspots OOS MSE, MAPE. QARIMA vs. classical. Fig. 6: Sunspots DM (MSE) p-values with α reference line; significant QARIMA bars appear in green. Fig. 7: Sunspots DM (MAE) − log 10 (p) with threshold at − log 10 (α). Interpretation of Sunpots results Within our pipeline quantum-informed differencing (Q-d estimation), quantum- assisted lag discovery for candidate (p,d,q), and VQC-only coefficient estimation under a fixed optimizer budget QARIMA retains the Box–Jenkins structure while replac- ing linear estimation with a shallow VQC (reps=1) that nonlinearly re-embeds lagged information before the ARIMA update. Even with a single repetition, the unitary embedding acts as a learned feature map over lags, expanding the hypothesis class beyond linear autoregression yet preserving parsimony in (p,d,q). For sunspots, a series with long-memory and quasi-periodic dynamics (the ∼11-year Schwabe cycle with drifting amplitude/phase), this shallow VQC, together with Q-d estimation, 32 Table 13: CO2 Classical Vs Quantum OOS CO2 Classical Vs Quantum OOS ModelNMSEMAPE Classical (pmdarima non-seasonal) (5, 1, 0)12078.372040.022805 Quantum (pdq=(10, 1, 1))12010.025750.007472 Quantum (pdq=(9, 1, 1))12010.176550.007496 Quantum (pdq=(8, 1, 1))12010.458980.007573 Quantum (pdq=(3, 1, 1))12079.37850.022943 Quantum (pdq=(4, 1, 1))12079.969090.023035 Quantum (pdq=(5, 1, 1))12081.124830.023214 Quantum (pdq=(6, 1, 1))12081.131410.023215 Quantum (pdq=(7, 1, 1))12083.661550.023586 captures higher-order lag interactions and slow amplitude/phase modulation with- out introducing an explicit seasonal operator. This explains why richer quantum orders (e.g., Q(10,1,1), Q(8/9,1,3)) generalize better across the 128-step OOS window, adapting to regime variation while keeping the interpretability provided by differenc- ing. Consequently, QARIMA yields statistically significant and practically meaningful reductions in OOS MSE/MAPE for most candidate models, consistent with Table 12 and the accompanying DM panels. 5.4 Muna Lua CO2 dataset We analyze the canonical Mauna Loa atmospheric CO 2 record (monthly ppm). The series features a strong upward trend (Keeling curve) and a pronounced annual cycle whose amplitude drifts slowly over time. Hence it is nonstationary in level with deterministic seasonality and gradual modulation, stressing purely linear, low–order ARIMA baselines. The Mauna Loa CO 2 series comprises monthly observations and after removing initial missing values it contains 468 observations in total. We reserve the last 120 months for out-of-sample (OOS) testing and use the preceding 348 months for training. The generated QARIMA models are evaluated against a classical non- seasonal baseline pmdarima ARIMA(5,1,0). OOS results for MSE and MAPE, together with DM tests (MSE/MAE), are reported in Table 13 and Figs. 8–10. Muna Lua CO2 Performance assessment (OOS). MSE MAPE comparison. Following the same protocol as Sunspots, we com- pare multiple QARIMA(p,d,q) models against a non–seasonal classical comparator ARIMA(5,1,0) on an OOS window of N =120 months. Table 13 and Fig. 8 show that high–p quantum models with a small MA term Q(10,1,1), Q(9,1,1), Q(8,1,1), achieve large gains over the classical baseline: MSE drops from ≈ 78.37 to ≈ 10 and MAPE from ≈ 2.28% to ≈ 0.75%. Lower–order quantum variants Q(3–7,1,1) do not improve upon the baseline (MSE ≈ 79–86, MAPE ≈ 2.29–2.36%), indicating that the CO 2 series benefits from richer autoregressive memory with modest q when using a non–seasonal specification. 33 Fig. 8: CO2 OOS MSE, MAPE. QARIMA vs. classical. Fig. 9: CO2 DM (MSE) p-values with α reference line; significant QARIMA bars appear in green. Fig. 10: CO2 DM (MAE) − log 10 (p) with threshold at − log 10 (α). Diebold–Mariano (DM) tests. DM results (Figs. 9–10) corroborate the error patterns. Where available, models with lower mean loss than the classical baseline (e.g., the strong Q(10/9/8,1,1) group) show significant differences under MSE/MAE DM (bars below α or above the − log 10 (α) threshold).Conversely, quantum models with higher mean loss (Q(3–7,1,1)) are also significantly different, but the direction (worse than classical) matches their positive classical–minus–quantum loss differ- entials. Complete per–model DM tables and mean–loss deltas are provided in the Appendix. 34 Table 14: AusBeer classical vs. quantum OOS (last 8 quarters). ModelNMSEMAPE Classical (pmdarima non-seasonal) (0, 1, 1)81491.7620630.080529 Quantum (pdq=(7, 1, 1))859.7923370.016178 Quantum (pdq=(10, 1, 1))876.1900980.017387 Quantum (pdq=(9, 1, 1))884.1308230.017418 Quantum (pdq=(6, 1, 1))895.9322820.018874 Quantum (pdq=(3, 1, 3))896.0296410.017407 Quantum (pdq=(5, 1, 1))899.1748030.019017 Quantum (pdq=(2, 1, 6))8143.4763360.024085 Quantum (pdq=(1, 1, 6))81848.8558560.074425 Interpretation of CO2 results. In line with our pipeline quantum-informed differencing (Q-d), quantum-assisted lag discovery for (p,d,q), and VQC-only coefficient estimation under a fixed opti- mizer budget with a shallow variational embedding (reps=1) QARIMA keeps the Box–Jenkins structure but replaces linear estimation by a unitary, learned feature map over lags. For CO 2 , this embedding enhances representation of long memory (via larger p) and introduces mild nonlinear mixing that effectively tracks the annual cycle’s amplitude drift without an explicit seasonal operator. As a result, high–p QARIMA models (Q(10/9/8,1,1)) deliver large, statistically significant OOS error reductions, while lower–order quantum variants (Q(3–7,1,1)) confirm via DM that the classical baseline is preferable when autoregressive memory is too limited. 5.5 Australian Beer Production Dataset The quarterly ausbeer series contains 211 observations. We follow the same protocol as for the other datasets and keep the last 8 quarters for out-of-sample (OOS) evaluation, using the first 203 quarters for training: N total = 211, N train = 203, N OOS = 8. All quantum ARIMA (QARIMA) models are compared against a simple classical non-seasonal baseline pmdarima ARIMA(0, 1, 1). 35 Fig. 11: AusBeer OOS MSE and MAPE (last 8 quarters). QARIMA vs. classical ARIMA(0,1,1). Fig. 12: AusBeer DM (MSE) p-values with α = 0.05 reference line; significant QARIMA models appear below the line. Fig. 13: AusBeer DM (MAE) − log 10 (p) with threshold at − log 10 (α). Australian beer OOS performance assessment. MSE / MAPE comparison. For the AusBeer series we used the last 8 quar- ters as the OOS window and compared all QARIMA candidates against the classical non-seasonal pmdarima ARIMA(0, 1, 1). The OOS results are given in Table 14 and visualised in Fig. 11. Unlike the short-history industrial series (Woolyarn), AusBeer shows a very strong separation: the classical ARIMA(0, 1, 1) records an OOS MSE 36 of 1491.8 and MAPE of 0.0805, whereas the best quantum model, Q(7, 1, 1), lowers these to 59.8 (MSE) and 0.0162 (MAPE) on the same 8-point horizon. A second tier of quantum models, Q(10, 1, 1), Q(9, 1, 1), Q(6, 1, 1), and Q(3, 1, 3) also stays clearly below 100 MSE and below 0.019 MAPE, showing that the improvement is shared by a family of quantum-configured orders and not a single outlier. Only the overparam- eterised MA-heavy quantum variants (e.g. Q(1, 1, 6)) drift back toward the classical error level. DM analysis. To confirm that these large numerical gains are not an artefact of the short 8-point window, we applied the Diebold–Mariano test against the classical ARIMA(0, 1, 1) forecast. In the MSE-based DM panels (Fig. 12) the leading quantum models, Q(7, 1, 1), Q(10, 1, 1), Q(9, 1, 1), Q(6, 1, 1), and Q(3, 1, 3) all lie well below the α = 0.05 line (and above the − log 10 (α) threshold), indicating that their forecast loss sequences are statistically different and, given their lower mean loss, better than the classical baseline. Models that are closer to the baseline (e.g. Q(2, 1, 6)) show weaker or no significance, which is consistent with their higher MSE/MAPE. On the MAE DM (Fig. 13) the pattern is similar: quantum models that already improved absolute error in Table 14 remain significant, while those that did not improve MAE fail to clear the α threshold. Together, the MSE/MAPE bars and DM panels make the case that, for AusBeer, quantum ARIMA’s lag search plus VQC estimation corrects a structural underfit in the classical ARIMA(0, 1, 1) baseline. Interpretation of AusBeer results. This outcome is fully consistent with the procedure described in Section X: we first set d using the quantum-informed differencing routine, then generate candidate lags from the quantum-inspired PACF/swap-test logic, and finally estimate the selected (p,d,q) models with a shallow VQC (reps=1). On AusBeer, this pipeline preferred higher AR orders (e.g. (7, 1, 1), (9, 1, 1), (10, 1, 1)), which are better aligned with the underlying quarterly pattern than the classical baseline ARIMA(0,1,1). Because the quantum layer re-embeds the selected lags before the ARIMA update, the model can simultaneously honor the differenced level and fit the short-horizon seasonal swing present in the last eight quarters. The very large and consistent error gap in Table 14 therefore reflects structural underfitting of the classical baseline, not instability of the DM test. 5.6 Australian Woolyran Production Dataset The woolyrnq series contains quarterly Australian woollen yarn production and is widely used as a benchmark for low–to–moderate volatility, industry–style time series. In our setting the series has N = 119 observations; we follow the same rolling/OOS protocol as in the previous datasets and keep the last 55 points for out-of-sample (OOS) evaluation, using the first 64 points as the history/training segment. This is therefore a short-history problem, in contrast to Sunspots and Mauna Loa CO 2 . All QARIMA models are evaluated against a classical non-seasonal pmdarima ARIMA(6, 1, 0) baseline. OOS results (MSE, MAPE) together with DM tests (MSE/MAE) are shown in Table 15 and Figs. 14–16. 37 Table 15: Woolyarn Classical vs Quantum OOS MSE & MAPE ModelNMSEMAPE Classical (pmdarima non-seasonal) (6, 1, 0)55528229.50640.105291 Quantum (pdq=(9, 1, 1))55533331.39270.104924 Quantum (pdq=(6, 1, 1))55530506.15270.105351 Quantum (pdq=(8, 1, 1))55538097.77940.105399 Quantum (pdq=(5, 1, 1))55544440.15890.107432 Quantum (pdq=(10, 1, 1))55555406.31540.107201 Quantum (pdq=(7, 1, 1))55564428.86580.108122 Quantum (pdq=(2, 1, 1))55579444.41060.112504 Quantum (pdq=(3, 1, 1))55594266.21210.113038 Quantum (pdq=(4, 1, 1))55609508.44870.113703 Quantum (pdq=(1, 1, 1))55554385.96390.112447 Fig. 14: Woolyarn OOS MSE (left) and MAPE (right): QARIMA vs. classical ARIMA(6,1,0). The classical model is best; several quantum models are close but not better. Woolyarn OOS Performance assessment. MSE/MAPE Comparision. Unlike Sunspots and CO 2 , the woollen yarn series does not reward deeper quantum orders. Fig. 14 shows that the classical ARIMA(6,1,0) achieves the lowest MSE (≈ 5.23× 10 5 ) and the lowest MAPE (≈ 0.1053). The closest quantum competitors are Q(6,1,1), Q(9,1,1), and Q(8,1,1), all of which stay in the narrow band 5.29–5.38×10 5 for MSE and 0.1049–0.1054 for MAPE. All other quantum variants (Q(5,1,1), Q(10,1,1), Q(7,1,1), Q(3,1,1), Q(4,1,1)) are clearly worse, with both MSE and MAPE climbing steadily as p is increased or the order is mismatched. This pattern is consistent with a short-history, low-noise, business–style series: the classical model already captures the available linear structure, and additional quantum flexibility can overfit the 64-point history. DM Test MSE/MAE. The DM–MSE panels (Fig. 15) show extremely small p- values for the high–order quantum models, e.g. Q(2,1,1), Q(1,1,1), Q(3,1,1), Q(4,1,1) reach− log 10 (p) in the range 10–15, far above the threshold− log 10 (0.05)≈ 1.30. This 38 Fig. 15: Woolyarn DM (MSE): almost all quantum models are significantly different from classical; the difference is mostly in the direction of higher loss. Fig. 16: Woolyarn DM (MAE): most quantum models remain significant; Q(6,1,1) and especially Q(5,1,1) are the only ones close to “no significant difference”. means their forecast loss sequences are clearly different from the classical baseline, and since their mean losses are larger, the inference is that classical ARIMA(6,1,0) is the more accurate model for this series. On the MAE DM (Fig. 16), most quantum models are again significant, but two models stand out: Q(6,1,1) has p ≈ 0.075 and Q(5,1,1) has p≈ 0.44, so for these two we cannot reject equal forecast accuracy under MAE at α = 0.05. This matches the OOS charts: these two are the only quantum models that sit close to the classical bar. Interpretation of Woolyarn results. This dataset is important because it shows the limits of the quantum approach under the very same pipeline (quantum-informed differencing for d, quantum-assisted lag discovery for (p,d,q), and shallow VQC-only coefficient estimation with reps=1). For a short history of 64 points, the linear structure learned by the classical ARIMA(6,1,0) appears to be already close to optimal; adding a unitary feature map over lags does not reveal substantial additional nonlinear or long-memory content. When we push p higher than what the series can support, the quantum models become distinguishable in DM (very small p) but mostly in the wrong direction, they are statistically different because they overfit or mis-specify. The only quantum models that behave competi- tively, Q(6,1,1) and Q(5,1,1), are precisely those that stay close to the classical order and add only a small MA component, which is in line with our manuscript’s design goal of capacity-controlled quantum ARIMA: keep the Box–Jenkins skeleton, let the 39 shallow VQC re-embed lags, but do not exceed what the data length can reliably esti- mate. Thus, the woollen yarn experiment complements the Sunspots and CO 2 results: QARIMA provides clear gains on long, nonlinear, or mixed-trend series, but on short, well-behaved industrial data the classical ARIMA can still be preferable, and the DM diagnostics make this explicit. 5.7 Sydney Weather 2024 dataset from NOAA 95768099999 For the weather use case, we extracted the Sydney 2024 summer segment (Decem- ber–February) from a NOAA Global Hourly station feed. The source was the NOAA NCEI endpoint [17], and we used the CSV file 95768099999.csv, which corresponds to a Sydney-area station in the specified bounding box. From this file we constructed a univariate daily temperature series and split it as history = 1782 observations and target/OOS = 336 observations, from which the 2024 summer slice (Dec–Feb) was taken. Table 16: Sydney 2024 Summer Temp: Classical vs. Quantum OOS (MSE & MAPE) ModelNMSEMAPE Classical (pmdarima non-seasonal) (2, 0, 1)33611.4359050.127427 Quantum (pdq=(3, 1, 1))33611.3562420.126633 Quantum (pdq=(4, 1, 1))33611.3934980.127034 Quantum (pdq=(5, 1, 1))33611.4120220.127216 Quantum (pdq=(6, 1, 1))33611.4234520.127318 Quantum (pdq=(9, 1, 1))33611.7026930.129946 Quantum (pdq=(8, 1, 1))33611.7359860.130267 Quantum (pdq=(10, 1, 1))33611.7430070.130273 Quantum (pdq=(7, 1, 1))33611.8936790.131706 Sydney weather temp OOS Performance assessment. MSE/MAPE Comparision. As with the other datasets, all quantum ARIMA (QARIMA) specifications were benchmarked against a classical, non-seasonal pmdarima baseline, which in this case selected ARIMA(2, 0, 1). Table 16 reports the OOS MSE and MAPE for this window. Over this summer slice the classical ARIMA(2, 0, 1) and most QARIMA runs achieve very similar errors. The shallow quan- tum variants Q(3,1,1), Q(4,1,1), Q(5,1,1), and Q(6,1,1) are marginally better (best: Q(3,1,1), MSE = 11.356, MAPE = 0.126653), but the improvement over the classi- cal model is modest. Higher-order quantum models Q(7,1,1), Q(8,1,1), Q(9,1,1), and Q(10,1,1) yield slightly larger MSE/MAPE yet remain in the same overall performance band, suggesting that on this short, low-variance horizon additional AR depth offers little incremental benefit rather than causing a clear deterioration (Fig. 17). 40 Fig. 17: Sydney Summer 2024 temperature: OOS MSE and MAPE, QARIMA vs. classical. Fig. 18: Sydney Summer 2024 temperature: DM (MSE) p-values and − log 10 (p) with α = 0.05. Q(4,1,1) and Q(3,1,1) are significant. DM analysis. To check whether the small numerical edges of the better quantum models are genuine, we applied the Diebold–Mariano (DM) test against the classi- cal ARIMA(2, 0, 1) forecast. On the MSE loss, the shallow quantum configurations Q(3,1,1) and Q(4,1,1) are the only ones that cross the α = 0.05 significance line in the p-value panels (see Figs. 18a–b), which is consistent with the fact that they are also the best (or near-best) models in the MSE/MAPE charts. The other quantum models, including the higher-AR variants (Q(7,1,1)–Q(10,1,1)), stay above the thresh- old: their errors are very close to the classical one and the short summer horizon does not provide enough observations for DM to declare a difference. For the MAE loss we additionally evaluated two MA-heavy quantum specifications, Q(1,1,6) and Q(2,1,6), which were included to test absolute-error behaviour. These two models yield clearly smaller p-values (DM–MAE p ≈ 5× 10 −4 and p ≈ 9× 10 −3 , respectively), and their − log 10 (p) bars rise well above the − log 10 (0.05) reference (Figs. 19a–b). All remain- ing models are non-significant on MAE, which matches the very narrow spread seen in the OOS MAPE chart. Interpretation of Sydney Summer Temperature Results. This seasonal weather slice illustrates that, on short and relatively low-variance hori- zons, our pipeline yields only modest but still detectable improvements. We followed the same procedure as in the manuscript consisting (i) quantum-informed selection of the differencing order, (i) quantum-assisted proposal of viable (p,d,q) candidates, and (i) shallow VQC-only estimation (reps=1) to embed the chosen lags before the 41 Fig. 19: Sydney Summer 2024 temperature: DM (MAE) p-values and− log 10 (p). MA- heavy Q(1,1,6) and Q(2,1,6) beat the classical baseline. ARIMA update. In this setting the quantum stage behaves mainly as a lag refiner : it steers the search toward compact models such as Q(3,1,1) and Q(4,1,1) that match the summer dynamics slightly better than the generic classical ARIMA(2, 0, 1), but the gain over the baseline is small, as the MSE/MAPE tables show. When MA-heavy quantum variants (Q(1,1,6), Q(2,1,6)) are included, the same pipeline can be aligned with an absolute-error objective, which explains why those two models stand out in the DM–MAE panels. Overall, the Sydney summer results show that QARIMA can register statistically visible effects even when the classical model is already strong and the evaluation window is short, and it does so without departing from the Box–Jenkins structure. 6 Discussion and Conclusion This manuscript introduced a quantum-informed procedure for univariate time-series forecasting that augments a classical ARIMA pipeline with lightweight quantum com- ponents for order discovery, lag refinement, and model selection. Instead of replacing the full statistical workflow, the proposed QARIMA framework targets the early steps that are often heuristic in practice choice of differencing, identification of viable (p,d,q) triples, and selection among models with similar likelihoods and makes them more systematic by leveraging quantum-inspired similarity tests and variational estimation. The resulting models remain interpretable in the classical ARIMA sense, require only shallow circuits, and can be trained with current simulators or small devices. Across five datasets of different character, Sunspots, Mauna Loa CO 2 , Australian woollen-yarn production, Australian beer production, and a 2024 Sydney weather slice, the quantum-augmented candidates consistently matched, and in several cases outperformed, the pmdarima-based baselines on out-of-sample MSE and/or MAPE, under the same train/OOS splits. The gains were clearest on series with slower or multi-scale dynamics (Mauna Loa, woollen yarn), where the quantum step tended to favour longer AR memories than the purely classical selection, and on short seasonal windows (Sydney 2024), where the quantum stage acted mainly as a lag refiner that nudged the solution toward compact Q(3,1,1)-Q(4,1,1) structures. On benchmark-like signals (Sunspots), improvements were smaller but still detectable when evaluated with Diebold-Mariano tests under MSE/MAE, confirming that the proposed procedure can differentiate between closely performing ARIMA variants on limited OOS blocks. 42 A practical outcome of the study is that the quantum part does not need to be deep or hardware-demanding to be useful: a single-layer VQC/swap-test stage, applied before the ARIMA fit, was enough to filter or re-rank candidate lag sets and to steer the model toward better OOS behaviour. This makes the approach deployable along- side existing Python time-series stacks (pmdarima, statsmodels) and compatible with hybrid classical-quantum experimentation. Future work will extend the same design to (i) richer seasonal and multivariate settings, (i) automated statistical testing pipelines so that DM significance is always reported with the metric, and (i) direct comparisons on NISQ hardware to validate the simulation-based results reported here. Acknowledgment The authors express gratitude to the IBM Quantum Experience platform and its team for creating the Qiskit platform and granting free access to their simulators for executing quantum circuits and conducting the experiments detailed below. The authors express appreciation for the Centre for Quantum Software and Information (CQSI) . 7 Statements and Declarations Competing Interests: The authors have no financial or non-financial competing interests. Authors’ contributions: The authors confirm their contribution to the paper as follows: Study conception and design: N.M., B.K.B., B.M., P.D.; Data collection: N.M.; Analysis and interpretation of results: N.M., B.K.B., B.M., P.D.; Draft manuscript preparation: N.M., B.K.B., B.M.; All authors reviewed the results and approved the final version of the manuscript. Funding: Authors declare that there has been no external funding. Availability of data and materials: All the data provided in this manuscript is generated during the simulation and can be provided upon reasonable request. References [1] Box GEP, Jenkins GM, Reinsel GC, Ljung GM. Time Series Analysis: Forecasting and Control. Wiley; 2015. [2] Akaike H. A new look at the statistical model identification. IEEE transactions on automatic control. 1974;19(6):716–723. [3] Schuld M, Banchi L, Gili G, Carleo G. Machine learning models that remember the past: quantum learning of non-Markovian processes. Physical Review A. 2021;104(2):022422. [4] Cerezo M, Arrasmith A, Babbush R, Benjamin SC, Endo S, Fujii K, et al. Variational quantum algorithms. Nature Reviews Physics. 2021;3:625–644. 43 [5] Lloyd S, Schuld M, Ijaz AW, Izaac J, Killoran N. Quantum embeddings for machine learning. arXiv preprint arXiv:200103622. 2020;. [6] Wilson GT. Review: Time Series Analysis: Forecasting and Control, 5th Edition by George E. P. Box, Gwilym M. Jenkins, Gregory C. Reinsel, and Greta M. Ljung. Journal of Time Series Analysis. 2016 Mar;37(5):709–711. https://doi. org/10.1111/jtsa.12194. [7] Buhrman H, Cleve R, Watrous J, de Wolf R. Quantum fingerprinting. Physical Review Letters. 2001;87(16):167902. [8] Schuld M, Killoran N. Quantum machine learning in feature Hilbert spaces. Physical Review Letters. 2019;122(4):040504. [9] Dickey DA, Fuller WA. Distribution of the Estimators for Autoregressive Time Series with a Unit Root. Journal of the American Statistical Association. 1979 Jun;74(366a):427–431. https://doi.org/10.1080/01621459.1979.10482531. [10] Shumway RH, Stoffer DS. Time Series Analysis and Its Applications: With R Examples. Springer; 2017. [11] Schuld M, Petruccione F.Supervised Learning with Quantum Computers. Quantum Science and Technology. Springer; 2018. [12] Powell MJD. A direct search optimization method that models the objective and constraint functions by linear interpolation. Advances in Optimization and Numerical Analysis. 1994;p. 51–67. https://doi.org/10.1007/978-94-015-8330-5 4. [13] Grant E, Wossnig L, Ostaszewski M, Benedetti M. An initialization strategy for addressing barren plateaus in parametrized quantum circuits. Quantum. 2019 Dec;3:214. https://doi.org/10.22331/q-2019-12-09-214. [14] Schuld M, Petruccione F. Machine Learning with Quantum Computers. Springer; 2021. [15] Havl ́ıˇcek V, C ́orcoles AD, Temme K, Harrow AW, Kandala A, Chow JM, et al. Supervised learning with quantum-enhanced feature spaces. Nature. 2019;567(7747):209–212. https://doi.org/10.1038/s41586-019-0980-2. [16] Diebold FX, Mariano RS. Comparing predictive accuracy. Journal of Business & Economic Statistics. 1995;13(3):253–263. [17] National Centers for Environmental Information (NCEI).: NOAA NCEI Global Hourly (ISD) — Search Results for DEW and WND within specified bound- ing box. Filtered search results (DEW, WND; bbox -33.839, 151.175, -33.903, 151.239; page 4).Available from: https://w.ncei.noaa.gov/access/search/ 44 data-search/global-hourly?dataTypes=DEW&dataTypes=WND&bbox=-33. 839,151.175,-33.903,151.239&pageNum=4. 45 8 Appendix 8.1 Sunspots Tables Table 17: Sunspots Classical VS Quantum - DM Stats MSE Sunspots Classical VS Quantum DM Stats MSE Quantum (pdq)blocksusedlossDMstatpvalueclassicalmeanlossquantummeanlossdeltameanloss (10, 1, 1)10MSE4.9417080.00000077442411.5317861563.023210848.508576 (9, 1, 3)10MSE4.8313550.00000135612411.5317861608.609316802.922470 (8, 1, 3)10MSE4.8120650.00000149382411.5317861613.363812798.167975 (6, 1, 3)10MSE4.6806260.00000286002411.5317861603.130658808.401129 (7, 1, 3)10MSE4.6084400.00000405702411.5317861629.559531781.972255 (5, 1, 3)10MSE4.4049790.00001057942411.5317861682.786868728.744918 (4, 1, 1)10MSE2.7547910.00587295502411.5317862056.829346354.702441 (3, 1, 1)10MSE1.9172760.05520290002411.5317862179.575810231.955977 Table 18: Sunspots Classical VS Quantum - DM Stats MAE Sunspots Classical VS QuantumDM Stats MAE Quantum pdqblocksusedlossDMstatpvalueclassicalmeanlossquantummeanlossdeltameanloss (10, 1, 1)10MAE4.9086690.000000917035.73411728.3876207.346497 (8, 1, 3)10MAE4.6092350.000004041535.73411729.1068196.627298 (9, 1, 3)10MAE4.6028180.000004168135.73411729.0665956.667522 (6, 1, 3)10MAE4.5393800.000005642035.73411728.8669356.867182 (7, 1, 3)10MAE4.3245370.000015285335.73411729.3069296.427188 (5, 1, 3)10MAE4.2426860.000022086035.73411729.8375745.896543 (4, 1, 1)10MAE1.1553930.247929500035.73411734.4239271.310189 (3, 1, 1)10MAE-0.0053230.995752600035.73411735.739900-0.005784 8.2 Woolyarn Tables Table 19: Woolyarn Classical VS Quantum - DM Stats MSE Woolyarn Classical VS QuantumDM Stats MSE Quantum pdqblocksusedlossDMstatpvalueclassicalmeanlossquantummeanlossdeltameanloss (2, 1, 1)1MSE-7.9683940.0000000000934889.8231136324-201434.0761 (1, 1, 1)1MSE-7.9520580.0000000000934889.8231185477-250586.7362 (3, 1, 1)1MSE-7.2454620.0000000000934889.8231131699-196809.3958 (4, 1, 1)1MSE-6.4276430.0000000001934889.8231042081-107190.8507 (7, 1, 1)1MSE-4.5000980.0000067922934889.8231343679-408789.1725 (10, 1, 1)1MSE-4.0711370.0000467843934889.8231026709-91818.85786 (8, 1, 1)1MSE-4.0236520.0000573025934889.8231024751-89861.61058 (9, 1, 1)1MSE-3.9371440.0000824572934889.8231009798-74908.25687 (6, 1, 1)1MSE-2.9924950.0027670760934889.823942247.1-7357.235698 (5, 1, 1)1MSE-2.5486270.0108148000934889.823948237.4-13347.5792 46 Table 20: Woolyarn Classical VS Quantum - DM Stats MAE Woolyarn Classical VS QuantumDM Stats MAE Quantum pdqblocksusedlossDMstatpvalueclassicalmeanlossquantummeanlossdeltameanloss (3, 1, 1)1MAE-7.369330.0000000000826.718879904.906818-78.187938 (4, 1, 1)1MAE-7.1485290.0000000000826.718879865.442915-38.724036 (2, 1, 1)1MAE-6.4341410.0000000001826.718879902.821988-76.103109 (1, 1, 1)1MAE-5.2205830.0000001784826.718879926.307162-99.588282 (7, 1, 1)1MAE-4.2725850.0000193220826.718879994.530279-167.811399 (10, 1, 1)1MAE-3.2510130.0011499450826.718879869.822437-43.103558 (8, 1, 1)1MAE-3.1944890.0014007860826.718879869.053863-42.334983 (9, 1, 1)1MAE-3.0395560.0023693500826.718879862.313859-35.594980 (6, 1, 1)1MAE-1.7788330.0752672000826.718879829.852377-3.133498 (5, 1, 1)1MAE0.7677600.4426297000826.718879823.7079483.010931 8.3 Muna Lua CO2 Tables Table 21: CO 2 Classical vs Quantum - DM Stats (MSE) CO 2 Classical vs Quantum - DM Stats (MSE) Quantum (pdq)blocksusedlossDMstatpvalueclassicalmeanlossquantummeanlossdeltameanloss (8, 1, 1)174MSE8.6811930.00E+009.0803823.5998585.480524 (9, 1, 1)172MSE8.9911400.00E+009.0993703.2827925.816578 (10, 1, 1)174MSE9.1942970.00E+009.0803823.1924905.887893 (5, 1, 1)174MSE-6.7299801.70E-119.0803829.236257-0.155875 (6, 1, 1)174MSE-6.6735772.50E-119.0803829.247513-0.167131 (4, 1, 1)174MSE-5.9479682.71E-099.0803829.217262-0.136880 (3, 1, 1)174MSE-5.2883341.23E-079.0803829.234463-0.154080 (7, 1, 1)174MSE5.2305771.69E-079.0803826.5636342.516748 Table 22: CO 2 Classical VS Quantum - DM Stats MAE CO 2 Classical VS QuantumDM Stats MAE Quantum pdqblocksusedlossDMstatpvalueclassicalmeanlossquantummeanlossdeltameanloss (8, 1, 1)174MAE9.2520650.00E+002.3157271.4060470.909681 (9, 1, 1)172MAE9.5613530.00E+002.3143381.3484580.965881 (10, 1, 1)174MAE9.8577190.00E+002.3157271.3248590.990868 (3, 1, 1)174MAE-5.6394721.71E-082.3157272.338754-0.023027 (7, 1, 1)174MAE5.4862844.10E-082.3157271.8877410.427986 (4, 1, 1)174MAE-5.4173996.05E-082.3157272.332440-0.016713 (5, 1, 1)174MAE-4.5363555.72E-062.3157272.328392-0.012664 (6, 1, 1)174MAE-4.4743347.66E-062.3157272.329375-0.013647 8.4 Australian Beer Production Tables 47 Table 23: AusBeer Classical VS Quantum - DM Stats (MSE) AusBeer Classical VS QuantumDM Stats MSE Quantum pdqblocksusedlossDMstatpvalueclassicalmeanlossquantummeanlossdeltameanloss (1, 1, 6)84MSE11.8780280.00E+002559.3285831296.1946171263.133966 (2, 1, 6)84MSE15.3530450.00E+002559.328583564.6066411994.721942 (7, 1, 1)84MSE15.6265420.00E+002559.328583388.1200112171.208572 (6, 1, 1)84MSE15.6269730.00E+002559.328583412.9911202146.337463 (5, 1, 1)84MSE15.7082020.00E+002559.328583408.2589032151.069680 (10, 1, 1)84MSE15.7306690.00E+002559.328583381.9987082177.329875 (9, 1, 1)84MSE15.7752090.00E+002559.328583378.5265932180.801990 (3, 1, 3)84MSE16.0626540.00E+002559.328583410.6654462148.663137 Table 24: AusBeer Classical VS Quantum - DM Stats (MAE) AusBeer Classical VS QuantumDM Stats MAE Quantum pdqblocksusedlossDMstatpvalueclassicalmeanlossquantummeanlossdeltameanloss (1, 1, 6)84MAE11.7723010.00E+0041.66672529.30605712.360668 (2, 1, 6)84MAE18.6321630.00E+0041.66672518.65199223.014733 (3, 1, 3)84MAE19.3755780.00E+0041.66672516.01142825.655297 (6, 1, 1)84MAE19.8592600.00E+0041.66672516.10459425.562131 (7, 1, 1)84MAE20.3027890.00E+0041.66672515.58837226.078353 (5, 1, 1)84MAE20.3274770.00E+0041.66672516.10933525.557390 (9, 1, 1)84MAE20.6688630.00E+0041.66672515.03750326.629222 (10, 1, 1)84MAE20.8544760.00E+0041.66672515.10646626.560259 8.5 Sydney weather 2024 tables Table 25: Sydney 2024 Summer Temp - Classical vs. QuantumDM Stats (MSE) Sydney 2024 Summer Temp Classical VS QuantumDM Stats MSE Quantum pdqblocksusedlossDMstatpvalueclassicalmeanlossquantummeanlossdeltameanloss (4, 1, 1)160MSE-3.419816.27E-048.16779.023377-0.855678 (3, 1, 1)161MSE-2.9434983.25E-038.17738319.809067-11.631684 (5, 1, 1)161MSE-2.4568011.40E-028.1773838.713587-0.536204 (6, 1, 1)161MSE-1.5448261.22E-018.1773838.668534-0.491151 (10, 1, 1)155MSE0.6106495.41E-018.2305467.9064640.324082 (7, 1, 1)161MSE-0.5774975.64E-018.1773838.416235-0.238851 (9, 1, 1)159MSE0.5553035.79E-018.2583267.9897770.268549 (8, 1, 1)161MSE0.1755058.61E-018.1773838.0996770.077707 Table 26: Sydney 2024 Summer Classical VS QuantumDM Stats (MAE) Sydney 2024 Summer Classical VS QuantumDM Stats MAE Quantum pdqblocksusedlossDMstatpvalueclassicalmeanlossquantummeanlossdeltameanloss (3, 1, 1)161MAE-4.0381515.40E-052.1837192.875055-0.691336 (4, 1, 1)160MAE-2.9758232.92E-032.1807982.280019-0.099221 (5, 1, 1)161MAE-1.7758337.58E-022.1837192.245571-0.061852 (10, 1, 1)155MAE1.2056822.28E-012.1934262.1063180.087108 (6, 1, 1)161MAE-1.1581402.47E-012.1837192.237213-0.053494 (9, 1, 1)159MAE1.0889502.76E-012.1982112.1244410.073770 (8, 1, 1)161MAE0.6019695.47E-012.1837192.1455490.038169 (7, 1, 1)161MAE-0.1640428.70E-012.1837192.193769-0.010050 48