arXiv CSopen access

QARIMA: A Quantum Approach To Classical Time Series Analysis

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

arXiv:2604.08277v1 [quant-ph] 9 Apr 2026

Nishikanta Mohanty1*, Bikash K. Behera2,3 , Badshah Mukherjee4 , Pravat Dash5 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 Università 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): [email protected]; Contributing authors: [email protected]; [email protected]; [email protected]; Abstract We present a quantum-inspired ARIMA methodology that integrates quantumassisted lag discovery with 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. 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 univariate sequences, since it decomposes signals into Auto-Regressive (AR), Differencing (I), and Moving-Average (MA) components. Although ARIMA models are interpretable and computationally efficient, they often fall short in capturing long-range dependencies, nonlinear dynamics, or stochastic effects that are inherent in realworld 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 highdimensional 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 similarity 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) differencing 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 regressors. 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 phasecorrected cosine similarities; candidate (p, q) are gated by these quantum diagnostics,

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 machinelearning techniques.

2 Background ARIMA models are among the most widely used approaches for modeling and forecasting 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 {yt }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: Bp (L) (1 − L)d yt = Θq (L) εt ,

(1)

where: • L is the lag operator, Lyt = yt−1 . • Bp (L) = 1 − b1 L − · · · − bp Lp denotes the AR polynomial. • Θq (L) = 1 + θ1 L + · · · + θq Lq 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: yt =

p X

bi yt−i + εt +

q X

θj εt−j .

(2)

j=1

i=1

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 crossvalidation. 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 optimization. These parameters will be further tuned by a VQC for improved performance in case-to-case scenarios.

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 projections with phase-corrected cosine similarities and feed directly into the delayed-matrix construction used for estimation. A delayed–matrix construction then aligns the quantum 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 outputs; 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 concrete 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 xt 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, θ ∈ Rn , we define normalized states:

∥x∥ ∥θ∥ ϕ = √ , −√ , Z Z

Z = ∥x∥2 + ∥θ∥2

1 x1 θ 1 xn θ n ψ=√ , ,..., , ∥x∥ ∥θ∥ 2 ∥x∥ ∥θ∥

(3)

(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 2p0 − 1,

p0 = Pr(measuring 0)

(5)

This estimation is executed using repeated quantum circuit simulations over a specified 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, θ ∈ Rn 1: Compute norms: ∥x∥, ∥θ∥ 2: if any norm is zero then 3: return default states [1, 0], [1, 0] 4: end if 2 2 5: Compute Z ← ∥x∥ + ∥θ∥ h i 6:

Compute ϕ ←

∥x∥ ∥θ∥ √ ,−√ Z Z

Initialize ψ ← [] for i = 1 to n do xi√ , θi√ 9: Append to ψ: ∥x∥ 2 ∥θ∥ 2 10: end for 11: return ϕ, ψ 7: 8:

Algorithm 2 Compact Swap Test Projection Require: Input x, θ ∈ Rn , number of shots S 1: Compute ϕ, ψ ← prep-swaptest(x, θ) 2: Normalize ϕ, ψ to unit ℓ2 norm k k 3: Pad ψ to length 2 for minimal k such that 2 ≥ 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 p0 √ 10: return cos(θ, x) ≈ 2p0 − 1

Algorithm 3 Compact Swap Test for Scaled Dot Product Require: Vectors x, θ ∈ Rn , 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: ϕ ← ϕ/∥ϕ∥, ψ ← ψ/∥ψ∥ ⌈log2 (len(ψ))⌉ 7: Compute target length: L ← 2 8: if len(ψ) < L then 9: Pad ψ with zeros to length L 10: end if 11: k ← log2 (L) ▷ Number of qubits for ψ 12: Allocate registers: one control qubit q1 , 1 + k ancilla qubits q2 , and 1 classical bit c 13: Initialize ϕ onto q2 [0] and ψ onto q2 [1], . . . , q2 [k] 14: Apply: X and H to control qubit q1 15: Apply: Controlled-swap between q2 [0] and q2 [1] 16: Apply: H to q1 and measure into c 17: Simulate circuit for S shots using QASM backend 18: Obtain measurement probabilities: p0 ← Pr(0), p1 ← Pr(1) 19: Compute swap test projection: Pswap ← 1 − 2p0 + p1 20:

return ∥x∥ · ∥θ∥ ·

max(Pswap , 0)

|qc ⟩

•

|ϕ⟩ CSWAP |ψ⟩ Fig. 1: Compact swap test estimating similarity between encoded states |ϕ⟩ and |ψ⟩. The control qubit is prepared in |0⟩, Hadamard → CSWAP → Hadamard, then measured to yield p0 , from which cos θswap and the entropy term H(p0 ) in LAR / LMA 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 traditional 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 generating differenced variants up to a pre-defined maximum order dmax (Algorithm 5). For each differencing level d, a pair of parameters (α, γ) are estimated to minimize the mean squared error (MSE) between the differenced target Deld and its predicted value ŷt , computed via swap-test–based quantum projection. Specifically, we define the projection-based prediction as ŷt = ⟨ψ(xt )|ψ(θ)⟩,

(6)

xt = [1, xt ] and θ = [α, γ]

(7)

where

are encoded into quantum states using amplitude encoding. The objective at each differencing level is to minimize: Ld (α, γ) =

1 X 2 (Deld (t) − ŷt ) . N t

(8)

The optimization is performed using a classical method such as COBYLA. A convergence 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.

Algorithm 4 Build Delay Matrix Require: Series y, maximum lag p 1: for i = 1 to p do 2: Create column lagi ← yt−i 3: end for 4: Append target yt return DataFrame of p lag columns and target yt

Algorithm 5 Generate Differenced Series Require: Series y, max differencing order dmax 1: Set prev ← y 2: for i = 1 to dmax do 3: Deli ← prev − prevt−1 4: prev ← Deli 5: end forreturn DataFrame with Del1 , . . . , Deldmax

Algorithm 6 Quantum-Inspired Estimation of Differencing Order d Require: Time series y, maximum differencing depth dmax , lag order p, convergence threshold ϵ, patience T , max optimizer iterations Imax , primary loss metric L, random seed s 1: Construct delay matrix from y using p lags (Algorithm 4) 2: Generate differenced series up to dmax (Algorithm 5) 3: Initialize d ← 0, metrics log ← ∅ 4: Evaluate d = 0: 5: Let xt ← lag1 (t), yt ← y(t) 6: Optimize (α0 , γ0 ) to minimize: L0 (α, γ) = 7:

1 X 2 (yt − ŷt ) N t

Compute predictions via swap test projection: ŷt = ⟨ψ([1, xt ]) | ψ([α, γ])⟩

Store α0 , γ0 , L0 in metrics log for d = 1 to dmax do 10: if Deld or Deld−1 not defined then 11: Break 12: end if 13: Let xt ← Deld−1 (t), yt ← Deld (t) 14: Initialize optimizer with seed s + d 15: Optimize (αd , γd ) to minimize: 8:

Ld (α, γ) = 16:

1 X 2 (yt − ŷt ) N t

Compute predictions via swap test projection: ŷt = ⟨ψ([1, xt ]) | ψ([α, γ])⟩

Store evaluation metrics in metrics log 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 Ld and γd ≈ 0 return d , αd∗ , γd∗

17:

18:

The estimation process is governed by several hyperparameters that control convergence behavior, optimization stability, and evaluation fidelity. These include the

maximum differencing depth dmax , maximum number of optimization iterations, convergence 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 Symbol

Name

Description

dmax p ϵ

Max differencing depth Lag order Convergence threshold

T Imax L

Patience Max optimizer iterations Loss metric

Random seed

Maximum number of differencing levels tested. Number of past lags used to construct the delay matrix. Threshold for detecting convergence of γ over recent iterations. Number of iterations required to confirm γ stability. Maximum steps allowed per optimization of (α, γ). Primary evaluation metric (e.g., MSE, MAE) used for model selection. Seed for reproducibility across differencing levels.

4.3 Quantum-Inspired Autocorrelation Function (ACF) The ACF quantifies linear dependence between yt and its lagged values yt−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 {yt }N t=1 (optionally centered and/or differenced), define aligned pairs {(yt , yt−k )}N for lag k. Let ψ(·) denote amplitude encoding of t=k+1 a real scalar into a single-qubit state (or a short register when batching); the swaptest estimate of cosine similarity between encodings is obtained from the control-qubit outcome probability p0 :

q [θt,k = cos

(t,k)

2p0

−1 ,

(9)

which serves as a quantum proxy for normalized correlation at lag k. Averaging over aligned pairs yields the quantum-inspired ACF value ρbkQ =

N X 1 [θt,k . cos N −k

(10)

t=k+1

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

(t,k) (t,k) (t,k) = cos θswap + ω θdot − θswap ,

(11)

P and set ρbkQ = N 1−k t cos(ϕcorr are then selected by thresholding |ρbkQ | t,k ). Significant lags√ using either (i) classical asymptotic bounds ±Z/ N (with Z = 1.96 for 95% bands), (ii) a percentile rule on magnitudes, or (iii) a mean+σ cutoff; a fallback band at βτ retains near-significant lags to stabilize downstream MA identification. In classical diagnostics, 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 y1:N (optionally centered/differenced), max lag K, shots S, threshold mode Tm ∈ {static, percentile, std}, phase weight ω, fallback ratio β Ensure: ACF values {ρbkQ }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 {(yt , yt−k )}N t=k+1 7: Initialize accumulator s ← 0, count m ← 0 8: for t = k + 1 to N do [θt,k 9: Encode yt , yt−k ; run compact swap test with S shots to get cos 10: Optionally compute classical cosine and apply phase correction to get cos(ϕcorr t,k ) [ 11: s ← s + cos(ϕcorr m←m+1 t,k ) or cos θ t,k ; 12: end for 13: ρbkQ ← s/m; Append ρbkQ to ACF 14: end for 15: Compute magnitudes ak ← ρ bkQ 16: if Tm = static then √ 17: τ ← Z/ N 18: else if Tm = percentile then 19: τ ← Percentile(a, q) 20: else if Tm = std then 21: τ ← mean(a) + σ std(a) 22: else 23: error 24: end if 25: τf ← β · τ ; Lsig ← {k : ak ≥ τ }; Lfb ← {k : τf ≤ ak < τ } 26: return ACF = {ρ bkQ }, L ← Lsig ∪ Lfb 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 Symbol

Name

Description

K S Tm Z q σ ω β

max lag shots threshold mode threshold z dynamic percentile std multiplier phase weight fallback ratio

Maximum lag evaluated for ACF. Number of swap-test measurements√per pair (yt , yt−k ). Thresholding strategy: static (±Z/ N ), percentile, or std. Z-score for static bounds (typically 1.96 for 95%). Percentile for magnitude cutoff (e.g., 90th). Multiplier for std-based bound (mean + σ·std). Weight for phase-corrected cosine blend (optional). Fallback fraction of the primary threshold (0 < β < 1).

natively with our swap-test losses used later in MA training, and (ii) 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 yt and yt−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 xt 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 =

N X 1 ⟨ψ(yt )|ψ(yt−k )⟩ , N −k

(12)

t=k+1

where ψ(·) denotes a quantum state encoding of the input scalar into qubit amplitude 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 based on classical asymptotic confidence intervals [10] N • 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

PACF while embedding quantum-aligned projection logic into lag selection, a technique 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 Tm , 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 yt ← {yk+1 , . . . , yN }, yt−k ← {y1 , . . . , yN −k } 7: for each t do 8: Compute projection: (k)

ρt

= ⟨ψ(yt )|ψ(yt−k )⟩ ≈ SwapTest(yt , yt−k , S)

end for P (k) Compute ρ̂k = N 1−k t ρt 11: Append ρ̂k to PACF 12: end for 13: Compute threshold τ : 14: if Tm = “static” then 15: τ ← √ZN 16: else if Tm = “percentile” then 17: τ ← top-q percentile of |PACF| 18: else if Tm = “std” then 19: τ ←µ+σ 20: end if 21: Define fallback threshold τf ← β · τ 22: Identify significant lags: 9:

10:

Ls = {k | |ρ̂k | ≥ τ }, 23:

Lf = {k | τf ≤ |ρ̂k | < τ }

Combine: Lfinal ← Ls ∪ Lf if f = True return PACF, Lfinal , τ

4.5 AutoRegressive Order p In ARIMA modeling, the autoregressive (AR) order p governs how many past values of a time series yt are used to predict its current value. Classical AR models rely on

Table 3: Hyperparameters for Quantum-Inspired PACF Estimation Symbol

Name

Description

K S Tm

max lag shots threshold mode

Z q

threshold numerator dynamic percentile

σ β f

std multiplier fallback ratio enable fallback

Maximum lag to compute PACF for. Number of measurement shots used in the swap test circuit. Strategy to determine the threshold: “static”, “percentile”, or “std”. Z-score for static thresholding (typically 1.96 for 95%). Percentile threshold (e.g., 90th percentile) for projection magnitude. Number of standard deviations above mean (used in “std” mode). Proportion of τ used as the fallback threshold. Boolean flag to include fallback lags in final selection.

linear regression, using dot-product formulations over lagged vectors. In our quantuminspired 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 process is formulated as follows. For a lag vector xt and coefficient vector b at time t, the classical prediction is: (dot)

ŷt

= x⊤ t b

(13)

We compute the classical cosine similarity as: cos θdot =

x⊤ t b̂ , ∥xt ∥

where b̂ =

b ∥b∥

(14)

The quantum-inspired variant estimates the cosine angle using a simulated swap test: cos θswap = SWAP(xt , 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: (quantum)

ŷt

= ∥xt ∥ · ∥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:

LAR =

(quantum)

yt − ŷt

+ λcos

X X (cos θdot − cos θswap )2 + λent H(cos θswap ) t

(17) 14

where the binary entropy is computed as:

H(p) = −p log2 p − (1 − p) log2 (1 − p)

(18)

The following algorithms compute the total loss LAR 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 {xt }, targets {yt }, candidate AR coefficients b, swap test function, hyperparameters λcos , λent , λphase Ensure: Total quantum-inspired AR loss LAR 1: Normalize b to obtain unit vector bunit 2: Pad bunit and xt to power-of-2 length 3: for each time step t do (dot) 4: Compute classical dot product: ŷt ← x⊤ t b 5: 6: 7: 8: 9: 10: 11:

x⊤ b

Compute classical cosine: cos θdot ← t∥xtunit ∥ Estimate cos θswap via compact swap test Compute angular correction: ∆θ ← θdot − θswap Corrected angle: θcorr ← θswap + λphase · ∆θ (quantum) ← ∥xt ∥ · ∥b∥ · cos(θcorr ) Quantum prediction: ŷt Entropy: Ht ← −p0 log2 p0 − (1 − p0 ) log2 (1 − p0 ) where p0 = 1 − cos2 θswap Accumulate: (quantum) 2

L += (yt − ŷt 12: 13:

) + λcos (cos θdot − cos θswap )2 + λent Ht

end for return L

Table 4: Hyperparameters for Quantum-Inspired AR Loss Estimation Symbol

Name

Description

shots

λcos

lambda cosine align

λent

lambda entropy

phase weight

noise std

padded dim

save csv

Number of swap test measurement shots used for cosine similarity estimation. Penalty weight applied to alignment difference between classical dot and swap-based cosine. Penalty weight for quantum entropy derived from swap test probability. Scaling factor to correct the phase discrepancy between classical and quantum cosine. Standard deviation of added Gaussian noise to simulate quantum fluctuations. Size of padded vector, adjusted to nearest power-of-two for swap test register. Boolean 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 LAR (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 retaining PACF interpretability. OLS initialization provides a warm start that reduces susceptibility to barren plateaus and poor local minima [13]. Let y ∈ RT denote the differenced series. For each p ∈ P, we form X(p) ∈ R(T −p)×p and solve: b(0) = arg min ∥y − X(p) b∥22 .

(19)

The initialized coefficients are embedded into the VQC as Ry (bj ) rotations. The VQC-based AR loss is defined as: LVQC AR (β) =

(quantum)

yt − ŷt

(β)

+ λcos

(cos θdot,t − cos θswap,t ) + λent

H(p0,t ) ,

(20) where the phase-corrected prediction is: (quantum)

ŷt

(β) = ∥xt ∥ ∥b(β)∥ cos(θswap,t + ω (θdot,t − θswap,t )) , 16

(21)

with:

b x⊤ t b(β) , ∥xt ∥ b(β) b , b(β) = ∥b(β)∥

cos θdot,t =

p0,t = 1 − cos2 θswap,t , H(p) = −p log2 p − (1 − p) log2 (1 − p).

(22) (23) (24) (25)

The optimal AR order is chosen as:

p⋆ = arg min LVQC b⋆(p) , AR p∈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 set P, number of VQC layers r, maximum iterations Tmax 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 LVQC AR (β) via COBYLA for Tmax iterations ⋆ 6: Store (p, b⋆(p) , LVQC AR (b(p) )) 7: end for VQC ⋆ ⋆ 8: p ← arg minp∈P LAR (b(p) ) ⋆ ⋆ 9: return p , b(p⋆ )

Table 5: Hyperparameters for VQC-Based AR Order Estimation Symbol

Name

Description

P d⋆ b(0) r S Tmax λcos , λent , ω

candidate lags best d b init vqc layers shots max iter loss weights

PACF-identified set of candidate AR orders. Optimal differencing order from quantum ADF tests. OLS-initialized AR coefficients. Number of VQC entangling layers. Number of measurement shots in swap tests. Maximum iterations for VQC optimization. Regularization and phase correction weights for the loss function.

(0)

Ry (β1 )

Ry b2

(0)

Ry (β2 )

.. .

Ry bp

Ry b1

(0)

(1)

•

(1)

(2)

Ry (β1 )

•

ŷt = ⟨O⟩

(2)

(1)

(2)

Ry (βp )

Fig. 2: VQC for AR order estimation with p qubits (one per lag). The first rotation (0) layer uses OLS initialization Ry (bi ); then r=2 variational layers (shown) alternate (ℓ) trainable Ry (βi ) with a linear entangling pattern (CNOT ladder). Readout yields ŷt = ⟨O⟩, and parameters are trained to minimize LAR (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 estimation, 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 = w1 , . . . , 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 bweak ∈ Rk denote the weak-lag coefficients and bfull = [ b⋆ ; bweak ] the concatenated vector. Using the quantum-inspired AR loss LAR (·) from Eq. (17), we optimize only bweak by minimizing Lweak (bweak ) = LAR bfull

+ λdev bweak − w 2 + λmag bweak 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

relative to a baseline loss Lbase (the best loss in the summary table) via λdev = 10−3 Lbase ,

λmag = 5 × 10−4 Lbase .

(28)

VQC parameterization: The weak-lag block is encoded by a depth-r VQC ansatz (e.g., Ry layers with entanglement). The rotation angles correspond to entries of bweak , while the anchor block b⋆ remains fixed. A classical derivative-free solver (COBYLA) refines bweak by minimizing Lweak , 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 Tmax Ensure: Extended order p′ , refined coefficients b⋆(p′ ) , refined loss Lweak ⋆ 1: Extract candidate magnitudes from Summary; remove magnitudes in |b |; keep k smallest as w = {w1 , . . . , wk } ′ ⋆ 2: Set p ← p + k ⋆ (p′ ) 3: Difference y by d and build AR design X and targets (0) (0) (0) ⋆ 4: Initialize bweak ← w and form bfull ← [ b ; bweak ] 5: Build a depth-r VQC with k trainable angles for bweak ⋆ 6: Define Lweak as in Eq. (27); fix b , optimize over bweak ⋆ 7: Minimize Lweak using COBYLA (budget Tmax ), obtain bweak ⋆ ⋆ ⋆ 8: Form b(p′ ) ← [ b ; bweak ] and evaluate Lweak ′ ⋆ 9: return p , b(p′ ) , Lweak

Table 6: Hyperparameters for VQC-Based Weak Lag Refinement Symbol

Name

Description

p⋆ b⋆ k w

best p anchor coeffs weak count weak inits

r Tmax λdev λmag

vqc depth max iter deviation penalty magnitude penalty

AR order chosen by PACF + VQC refinement. Fixed anchor coefficients for AR(p⋆ ). Number of weak lags added to form p′ = p⋆ + k. Weak-lag initialization magnitudes (from summary, excluding |b⋆ |). Number of entangling layers in the weak-lag VQC block. Max COBYLA iterations for minimizing Lweak . Weight on ∥bweak − w∥22 (scaled by Lbase ). Weight on ∥bweak ∥1 (scaled by Lbase ).

Anchor block (frozen)

(1)

Ry (βw,1 ) Weak-lag block (trainable)

(2)

•

Ry (βw,1 )

(1)

Ry (βw,2 )

(1)

Ry (βw,3 )

•

(2)

Ry (βw,2 )

(2)

Ry (βw,3 )

Fig. 3: VQC for weak-lag refinement with p⋆ frozen anchors (top block) and k trainable (p⋆ ) weak-lag qubits (bottom block). Anchors are initialized as Ry bi and kept fixed (ℓ) (shown as a grouped “frozen” block). Only weak-lag qubits carry trainable Ry (βw,j ) and are entangled (CNOT ladder restricted to the weak block). Training minimizes Lweak = LAR + λdev ∥bweak − w∥22 + λmag ∥bweak ∥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 = {w1 , w2 , . . . , wk } 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 h ⋆ i (p) binit = b(p ) , w1 , . . . , wj . We employ a VQC-based refinement over the weak-lag block, keeping the anchor coefficients fixed. The weak-lag optimization problem is formulated as:

2 (p) (p) (p) Lweak = LAR b(p) + λdev bweak − w1:j + λmag bweak 2

where: 20

, 1

(29)

• LAR 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 Lbase = LAR (b(p ) ): λdev = 10−3 Lbase ,

λmag = 5 × 10−4 Lbase .

(30)

For each refinement step, the weak-lag parameters are embedded into a j-qubit VQC with Ry rotations and r entangling layers, and optimized using COBYLA. The process continues progressively until the loss improvement (p)

∆(p) = Lweak − Lbase 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 Tmax ; threshold τ Ensure: Sequence {(p, b(p) , L(p) )} (p⋆ ) 1: Lbase ← LAR (b ) 2: for j = 1 to |w| do 3: p ← p⋆ + j ⋆ (p) 4: Initialize binit ← [b(p ) , w1 , . . . , wj ] 5: Build a j-qubit VQC (depth r) for the weak-lag block 6: Minimize L(p) (Eq. (29)) via COBYLA for Tmax iterations (p) 7: Record (p, b(p) , Lweak ) (p) 8: if Lweak − Lbase > τ then 9: break 10: end if 11: end for

Table 7: Hyperparameters for Progressive VQC-Based Weak Lag Refinement Symbol

Name

Description

w r Tmax λdev , λmag

weak lags vqc depth max iter penalty weights

Ordered list of weak-lag initialization values. Number of entangling layers in weak-lag VQC. Maximum COBYLA iterations per refinement step. Adaptive regularization weights based on Lbase .

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 residuals. These residuals represent the error between predicted and observed values of the differenced time series ∆d yt and are essential for both evaluating AR model accuracy and informing subsequent MA model design. For each lag order p, a delay matrix Xt is constructed from p prior values of the differenced series ∆d yt . The quantum-inspired prediction for each time step is formulated as:

ŷt = ∥b∥ · ∥Xt ∥ · cos (θt ) Xt · b̂ ∥Xt ∥ b b̂ = ∥b∥

cos(θt ) =

(31) (32) (33)

The residuals are then computed as: εt = ∆d yt − 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.

Algorithm 14 Residual Evaluation for Quantum-Inspired AR Models Require: Differenced time series ∆d y, optimal differencing order d, set of AR orders {pi } and coefficients {bi } Ensure: Residual statistics for each AR model 1: Initialize residual log R ← [ ] 2: for each row i in AR results do 3: p ← pi , b ← bi 4: Construct delay matrix X using p lags from ∆d y 5: Set target vector ytarget ← ∆d yp: 6: Normalize b to get unit vector b̂ 7: for each row xt in X do 8: Compute norm ∥xt ∥ xt ·b̂ 9: Compute cosine projection: cos(θt ) ← ∥x t∥ 10: Compute prediction: ŷt = ∥b∥ · ∥xt ∥ · cos(θt ) 11: Compute residual: εt = yt − 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 {yt }Tt=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 θ ∈ Rq 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.

4.7.1 Quantum MA Loss Estimation To estimate the quality of MA coefficients θ, we define a quantum-regularized loss function LMA that integrates classical residual prediction with quantum projections derived from swap tests [14, 15]. This formulation enables us to model the noise correlation 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)

ŷt = ∥θ∥ · ∥ϵt ∥ · cos(ϕcorrected ) t

(36)

and the predicted target as:

Here, ϕcorrected is a phase-adjusted angle between the normalized θ and ϵt , combining t classical and quantum projections: cos(ϕcorrected ) = cos (θswap + ω(θdot − θswap )) t

(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: LMA =

(yt − ŷt ) + λcos

X 2 Ht + λL2 ∥θ∥2 cos(ϕt ) − cos(ϕcorrected ) + λent t t

(38) where Ht is the entropy of the swap-test-derived distribution: Ht = −

pi log2 (pi ),

), with p0 = cos2 (ϕswap t

p1 = 1 − p0

(39)

Algorithm 15 Quantum-Inspired MA Loss Estimation Require: MA coefficients θ, residual matrix E, target series {yt }, trial ID, regularization weights λcos , λent , λL2 , phase weight ω, number of quantum shots S Ensure: Quantum-regularized loss LMA 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 Ht 10: Accumulate total loss components 11: end for 2 12: Add λL2 · ∥θ∥ regularization 13: return LMA

Table 8: Hyperparameters for Quantum MA Loss Symbol

Name

Description

λcos

lambda cosine align

λent

lambda entropy

ω S λL2

phase weight shots lambda l2

Penalty for misalignment between classical dot-product and quantum cosine. Weight for entropy regularization based on swap test measurements. Phase correction factor for angular adjustment. Number of measurements used in the swap test simulation. ℓ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 LMA across a candidate set Q = {1, 2, . . . , qmax }. 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)

q ⋆ = arg min LMA q∈Q

(40)

Each candidate q generates a delay matrix from residuals, against which the prediction 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, Iq )

(41)

The loss function LMA 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 {yt }, candidate range Q = [qmin , qmax ], loss function LMA Ensure: Optimal MA order q ⋆ 1: Initialize empty summary list 2: for each q ∈ Q do 3: Construct delay matrix Eq using q lags 4: Initialize MA vector θ (q) ∼ N (0, Iq ) (q) 5: Compute loss: LMA ← QuantumMAProjectionLoss(θ (q) , Eq , {yt }) (q) 6: Store (q, LMA ) in summary 7: end for (q) ⋆ 8: Select q ← arg min LMA ⋆ 9: return q

Table 9: Hyperparameters for MA Order Estimation Symbol

Name

Description

qmin , qmax θ init LMA

q range theta init loss function

Range of candidate MA orders tested during estimation. Random initialization of coefficient vector for each trial q. Quantum-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 LMA (θ) with 26

⋆

respect to the coefficient vector θ ∈ Rq . 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 LMA (θ) θ

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 build vqc ma(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 Ry rotation angles. The optimization 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 yMA . We then solve ⋆

θ (0) = arg min yMA − E(q ) θ 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 standard deviation of the residuals to quantify the remaining volatility in the MA model: v u T u1 X 2 (yt − ŷt ) , (44) σMA = t T t=1 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.

Algorithm 17 VQC-Refined MA Training for q ⋆ Require: Residuals {ϵt }, AR target {yt }, optimal order q ⋆ , VQC depth r Ensure: Optimized coefficients θ ⋆ , MA noise standard deviation σMA (q ⋆ ) 1: Construct delayed error matrix E from {ϵt } (0) 2: Initialize θ via OLS; fallback to U (−1, 1) ⋆ 3: Embed θ into build vqc ma(q , r) 4: Minimize LMA (θ) via COBYLA with ∥θ∥2 ≤ τ q P 1 2 5: Compute ϵt = yt − ŷt and σMA = t ϵt T ⋆ 6: return θ , σMA

Table 10: Hyperparameters for VQC-Refined MA Training Symbol (0)

θ τ r Tmax S λcos , λent , λL2 , ω

Name

Description

theta init norm constraint vqc reps max iter shots loss weights

Initial coefficient vector from OLS (clipped) or U (−1, 1). Norm bound on θ to control overfitting. Number of entangling layers in the VQC ansatz. Maximum iterations for COBYLA. Number of swap-test shots for loss evaluation. Loss weights and phase-correction factor (as in Eq. 38).

Ry (θ0 )

Ry (θ1 )

Ry (θ2 )

···

• •

··· ···

.. . qq⋆ −2

Ry (θq⋆ −2 )

qq⋆ −1

Ry (θq⋆ −1 )

•

··· ···

Fig. 4: Schematic VQC ansatz for MA coefficient refinement build vqc ma(q ⋆ , r). Each qubit qi receives a single-qubit Ry (θ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 LMA .

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 Column

Type

Description

p d q b θ σAR σMA

Integer Integer Integer Vector Vector Float Float

AR order Differencing order MA order (quantum-optimized) AR coefficients MA coefficients Noise std from AR residuals Estimated 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 CO2 , Australian Woollen Yarn Production, Australian Beer Production, and Sydney 2024 Weather. For each dataset, we followed the Methods pipeline: (i) identical preprocessing and differencing, (ii) generation of candidate QARIMA(p, d, q) models via the manuscript procedure, and (iii) 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).

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 provides 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 = MSEclassical − MSEquantum ,

∆MAPE = MAPEclassical − MAPEquantum .

Individual dataset subsections interpret these numbers in context (trend-like vs. seasonal 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 pvalue, and mark results that satisfy p ≤ α (typically α = 0.05). Alongside significance, we also show the corresponding classical and quantum mean losses and their difference ∆ = lossclassical − lossquantum , 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 (ii) achieves DM significance (p ≤ 0.05) under at least one loss (preferably both). When multiple QARIMA specifications 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 We use the classic annual mean sunspot counts provided in 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 modulation) 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 Model N M SE Classical (pmdarima) (2, 0, 0) seasonal=None 128 2181.589 Quantum (pdq=(3, 1, 1)) 128 2146.926 Quantum (pdq=(10, 1, 1)) 128 2160.998 Quantum (pdq=(7, 1, 3)) 128 2174.622 Quantum (pdq=(4, 1, 1)) 128 2179.405 Quantum (pdq=(9, 1, 3)) 128 2182.978 Quantum (pdq=(8, 1, 3)) 128 2184.006 Quantum (pdq=(6, 1, 3)) 128 2210.434 Quantum (pdq=(5, 1, 3)) 128 2256.788

MAPE 1.661 1.790 1.534 1.629 1.705 1.601 1.580 1.349 1.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 specifications 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 relative 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 variability, 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 occasionally 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) − log10 (p) with threshold at − log10 (α). Interpretation of Sunpots results Within our pipeline quantum-informed differencing (Q-d estimation), quantumassisted lag discovery for candidate (p, d, q), and VQC-only coefficient estimation under a fixed optimizer budget QARIMA retains the Box–Jenkins structure while replacing 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,

Table 13: CO2 Classical Vs Quantum OOS CO2 Classical Vs Quantum OOS Model Classical (pmdarima non-seasonal) (5, 1, 0) Quantum (pdq=(10, 1, 1)) Quantum (pdq=(9, 1, 1)) Quantum (pdq=(8, 1, 1)) Quantum (pdq=(3, 1, 1)) Quantum (pdq=(4, 1, 1)) Quantum (pdq=(5, 1, 1)) Quantum (pdq=(6, 1, 1)) Quantum (pdq=(7, 1, 1))

N 120 120 120 120 120 120 120 120 120

MSE 78.37204 10.02575 10.17655 10.45898 79.3785 79.96909 81.12483 81.13141 83.66155

MAPE 0.022805 0.007472 0.007496 0.007573 0.022943 0.023035 0.023214 0.023215 0.023586

captures higher-order lag interactions and slow amplitude/phase modulation without 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 differencing. 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 CO2 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 CO2 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 nonseasonal 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 compare 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 CO2 series benefits from richer autoregressive memory with modest q when using a non–seasonal specification.

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) − log10 (p) with threshold at − log10 (α). 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 − log10 (α) 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 differentials. Complete per–model DM tables and mean–loss deltas are provided in the Appendix.

Table 14: AusBeer classical vs. quantum OOS (last 8 quarters).

Model Classical (pmdarima non-seasonal) (0, 1, 1) Quantum (pdq=(7, 1, 1)) Quantum (pdq=(10, 1, 1)) Quantum (pdq=(9, 1, 1)) Quantum (pdq=(6, 1, 1)) Quantum (pdq=(3, 1, 3)) Quantum (pdq=(5, 1, 1)) Quantum (pdq=(2, 1, 6)) Quantum (pdq=(1, 1, 6))

N 8 8 8 8 8 8 8 8 8

MSE 1491.762063 59.792337 76.190098 84.130823 95.932282 96.029641 99.174803 143.476336 1848.855856

MAPE 0.080529 0.016178 0.017387 0.017418 0.018874 0.017407 0.019017 0.024085 0.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 optimizer 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 CO2 , 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: Ntotal = 211,

Ntrain = 203,

NOOS = 8.

All quantum ARIMA (QARIMA) models are compared against a simple classical non-seasonal baseline pmdarima ARIMA(0, 1, 1).

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) − log10 (p) with threshold at − log10 (α).

Australian beer OOS performance assessment. MSE / MAPE comparison. For the AusBeer series we used the last 8 quarters 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

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 overparameterised 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 − log10 (α) 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 CO2 . 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.

Table 15: Woolyarn Classical vs Quantum OOS MSE & MAPE

Model Classical (pmdarima non-seasonal) (6, 1, 0) Quantum (pdq=(9, 1, 1)) Quantum (pdq=(6, 1, 1)) Quantum (pdq=(8, 1, 1)) Quantum (pdq=(5, 1, 1)) Quantum (pdq=(10, 1, 1)) Quantum (pdq=(7, 1, 1)) Quantum (pdq=(2, 1, 1)) Quantum (pdq=(3, 1, 1)) Quantum (pdq=(4, 1, 1)) Quantum (pdq=(1, 1, 1))

N 55 55 55 55 55 55 55 55 55 55 55

MSE 528229.5064 533331.3927 530506.1527 538097.7794 544440.1589 555406.3154 564428.8658 579444.4106 594266.2121 609508.4487 554385.9639

MAPE 0.105291 0.104924 0.105351 0.105399 0.107432 0.107201 0.108122 0.112504 0.113038 0.113703 0.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 CO2 , 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 × 105 ) 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×105 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 pvalues 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 − log10 (p) in the range 10–15, far above the threshold − log10 (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 competitively, 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

shallow VQC re-embed lags, but do not exceed what the data length can reliably estimate. Thus, the woollen yarn experiment complements the Sunspots and CO2 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 (December–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)

Model Classical (pmdarima non-seasonal) (2, 0, 1) Quantum (pdq=(3, 1, 1)) Quantum (pdq=(4, 1, 1)) Quantum (pdq=(5, 1, 1)) Quantum (pdq=(6, 1, 1)) Quantum (pdq=(9, 1, 1)) Quantum (pdq=(8, 1, 1)) Quantum (pdq=(10, 1, 1)) Quantum (pdq=(7, 1, 1))

N 336 336 336 336 336 336 336 336 336

MSE 11.435905 11.356242 11.393498 11.412022 11.423452 11.702693 11.735986 11.743007 11.893679

MAPE 0.127427 0.126633 0.127034 0.127216 0.127318 0.129946 0.130267 0.130273 0.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 quantum 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 classical 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).

Fig. 17: Sydney Summer 2024 temperature: OOS MSE and MAPE, QARIMA vs. classical.

Fig. 18: Sydney Summer 2024 temperature: DM (MSE) p-values and − log10 (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 classical 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 threshold: 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 − log10 (p) bars rise well above the − log10 (0.05) reference (Figs. 19a–b). All remaining 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 horizons, 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, (ii) quantum-assisted proposal of viable (p, d, q) candidates, and (iii) shallow VQC-only estimation (reps=1) to embed the chosen lags before the

Fig. 19: Sydney Summer 2024 temperature: DM (MAE) p-values and − log10 (p). MAheavy 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 components 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 CO2 , 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.

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 alongside 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, (ii) automated statistical testing pipelines so that DM significance is always reported with the metric, and (iii) 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ı́ček V, Córcoles 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 bounding box. Filtered search results (DEW, WND; bbox -33.839, 151.175, -33.903, 151.239; page 4). Available from: https://www.ncei.noaa.gov/access/search/

data-search/global-hourly?dataTypes=DEW&dataTypes=WND&bbox=-33. 839,151.175,-33.903,151.239&pageNum=4.

8 Appendix 8.1 Sunspots Tables

Table 17: Sunspots Classical VS Quantum - DM Stats MSE Quantum (pdq) (10, 1, 1) (9, 1, 3) (8, 1, 3) (6, 1, 3) (7, 1, 3) (5, 1, 3) (4, 1, 1) (3, 1, 1)

blocks used 10 10 10 10 10 10 10 10

loss MSE MSE MSE MSE MSE MSE MSE MSE

Sunspots Classical VS Quantum DM Stats MSE DM stat p value classical mean loss quantum mean loss 4.941708 0.0000007744 2411.531786 1563.023210 4.831355 0.0000013561 2411.531786 1608.609316 4.812065 0.0000014938 2411.531786 1613.363812 4.680626 0.0000028600 2411.531786 1603.130658 4.608440 0.0000040570 2411.531786 1629.559531 4.404979 0.0000105794 2411.531786 1682.786868 2.754791 0.0058729550 2411.531786 2056.829346 1.917276 0.0552029000 2411.531786 2179.575810

delta mean loss 848.508576 802.922470 798.167975 808.401129 781.972255 728.744918 354.702441 231.955977

Table 18: Sunspots Classical VS Quantum - DM Stats MAE Quantum pdq (10, 1, 1) (8, 1, 3) (9, 1, 3) (6, 1, 3) (7, 1, 3) (5, 1, 3) (4, 1, 1) (3, 1, 1)

blocks used 10 10 10 10 10 10 10 10

Sunspots Classical VS Quantum DM Stats MAE loss DM stat p value classical mean loss quantum mean loss MAE 4.908669 0.0000009170 35.734117 28.387620 MAE 4.609235 0.0000040415 35.734117 29.106819 MAE 4.602818 0.0000041681 35.734117 29.066595 MAE 4.539380 0.0000056420 35.734117 28.866935 MAE 4.324537 0.0000152853 35.734117 29.306929 MAE 4.242686 0.0000220860 35.734117 29.837574 MAE 1.155393 0.2479295000 35.734117 34.423927 MAE -0.005323 0.9957526000 35.734117 35.739900

delta mean loss 7.346497 6.627298 6.667522 6.867182 6.427188 5.896543 1.310189 -0.005784

8.2 Woolyarn Tables

Table 19: Woolyarn Classical VS Quantum - DM Stats MSE Quantum pdq (2, 1, 1) (1, 1, 1) (3, 1, 1) (4, 1, 1) (7, 1, 1) (10, 1, 1) (8, 1, 1) (9, 1, 1) (6, 1, 1) (5, 1, 1)

blocks used 1 1 1 1 1 1 1 1 1 1

Woolyarn Classical VS Quantum DM Stats MSE loss DM stat p value classical mean loss quantum mean loss MSE -7.968394 0.0000000000 934889.823 1136324 MSE -7.952058 0.0000000000 934889.823 1185477 MSE -7.245462 0.0000000000 934889.823 1131699 MSE -6.427643 0.0000000001 934889.823 1042081 MSE -4.500098 0.0000067922 934889.823 1343679 MSE -4.071137 0.0000467843 934889.823 1026709 MSE -4.023652 0.0000573025 934889.823 1024751 MSE -3.937144 0.0000824572 934889.823 1009798 MSE -2.992495 0.0027670760 934889.823 942247.1 MSE -2.548627 0.0108148000 934889.823 948237.4

delta mean loss -201434.0761 -250586.7362 -196809.3958 -107190.8507 -408789.1725 -91818.85786 -89861.61058 -74908.25687 -7357.235698 -13347.5792

Table 20: Woolyarn Classical VS Quantum - DM Stats MAE Quantum pdq (3, 1, 1) (4, 1, 1) (2, 1, 1) (1, 1, 1) (7, 1, 1) (10, 1, 1) (8, 1, 1) (9, 1, 1) (6, 1, 1) (5, 1, 1)

blocks used 1 1 1 1 1 1 1 1 1 1

Woolyarn Classical VS Quantum DM Stats MAE loss DM stat p value classical mean loss quantum mean loss MAE -7.36933 0.0000000000 826.718879 904.906818 MAE -7.148529 0.0000000000 826.718879 865.442915 MAE -6.434141 0.0000000001 826.718879 902.821988 MAE -5.220583 0.0000001784 826.718879 926.307162 MAE -4.272585 0.0000193220 826.718879 994.530279 MAE -3.251013 0.0011499450 826.718879 869.822437 MAE -3.194489 0.0014007860 826.718879 869.053863 MAE -3.039556 0.0023693500 826.718879 862.313859 MAE -1.778833 0.0752672000 826.718879 829.852377 MAE 0.767760 0.4426297000 826.718879 823.707948

delta mean loss -78.187938 -38.724036 -76.103109 -99.588282 -167.811399 -43.103558 -42.334983 -35.594980 -3.133498 3.010931

8.3 Muna Lua CO2 Tables Table 21: CO2 Classical vs Quantum - DM Stats (MSE) CO2 Classical vs Quantum - DM Stats (MSE) Quantum (pdq) blocks used loss DM stat (8, 1, 1) 174 MSE 8.681193 (9, 1, 1) 172 MSE 8.991140 (10, 1, 1) 174 MSE 9.194297 (5, 1, 1) 174 MSE -6.729980 (6, 1, 1) 174 MSE -6.673577 (4, 1, 1) 174 MSE -5.947968 (3, 1, 1) 174 MSE -5.288334 (7, 1, 1) 174 MSE 5.230577

p value 0.00E+00 0.00E+00 0.00E+00 1.70E-11 2.50E-11 2.71E-09 1.23E-07 1.69E-07

classical mean loss 9.080382 9.099370 9.080382 9.080382 9.080382 9.080382 9.080382 9.080382

quantum mean loss 3.599858 3.282792 3.192490 9.236257 9.247513 9.217262 9.234463 6.563634

delta mean loss 5.480524 5.816578 5.887893 -0.155875 -0.167131 -0.136880 -0.154080 2.516748

Table 22: CO2 Classical VS Quantum - DM Stats MAE Quantum pdq (8, 1, 1) (9, 1, 1) (10, 1, 1) (3, 1, 1) (7, 1, 1) (4, 1, 1) (5, 1, 1) (6, 1, 1)

blocks used 174 172 174 174 174 174 174 174

loss MAE MAE MAE MAE MAE MAE MAE MAE

CO2 Classical VS Quantum DM Stats MAE DM stat p value classical mean loss quantum mean loss 9.252065 0.00E+00 2.315727 1.406047 9.561353 0.00E+00 2.314338 1.348458 9.857719 0.00E+00 2.315727 1.324859 -5.639472 1.71E-08 2.315727 2.338754 5.486284 4.10E-08 2.315727 1.887741 -5.417399 6.05E-08 2.315727 2.332440 -4.536355 5.72E-06 2.315727 2.328392 -4.474334 7.66E-06 2.315727 2.329375

8.4 Australian Beer Production Tables

delta mean loss 0.909681 0.965881 0.990868 -0.023027 0.427986 -0.016713 -0.012664 -0.013647

Table 23: AusBeer Classical VS Quantum - DM Stats (MSE) Quantum pdq (1, 1, 6) (2, 1, 6) (7, 1, 1) (6, 1, 1) (5, 1, 1) (10, 1, 1) (9, 1, 1) (3, 1, 3)

blocks used 84 84 84 84 84 84 84 84

AusBeer Classical VS Quantum DM Stats MSE loss DM stat p value classical mean loss quantum mean loss MSE 11.878028 0.00E+00 2559.328583 1296.194617 MSE 15.353045 0.00E+00 2559.328583 564.606641 MSE 15.626542 0.00E+00 2559.328583 388.120011 MSE 15.626973 0.00E+00 2559.328583 412.991120 MSE 15.708202 0.00E+00 2559.328583 408.258903 MSE 15.730669 0.00E+00 2559.328583 381.998708 MSE 15.775209 0.00E+00 2559.328583 378.526593 MSE 16.062654 0.00E+00 2559.328583 410.665446

delta mean loss 1263.133966 1994.721942 2171.208572 2146.337463 2151.069680 2177.329875 2180.801990 2148.663137

Table 24: AusBeer Classical VS Quantum - DM Stats (MAE) Quantum pdq (1, 1, 6) (2, 1, 6) (3, 1, 3) (6, 1, 1) (7, 1, 1) (5, 1, 1) (9, 1, 1) (10, 1, 1)

blocks used 84 84 84 84 84 84 84 84

AusBeer Classical VS Quantum DM Stats MAE loss DM stat p value classical mean loss quantum mean loss MAE 11.772301 0.00E+00 41.666725 29.306057 MAE 18.632163 0.00E+00 41.666725 18.651992 MAE 19.375578 0.00E+00 41.666725 16.011428 MAE 19.859260 0.00E+00 41.666725 16.104594 MAE 20.302789 0.00E+00 41.666725 15.588372 MAE 20.327477 0.00E+00 41.666725 16.109335 MAE 20.668863 0.00E+00 41.666725 15.037503 MAE 20.854476 0.00E+00 41.666725 15.106466

delta mean loss 12.360668 23.014733 25.655297 25.562131 26.078353 25.557390 26.629222 26.560259

8.5 Sydney weather 2024 tables

Table 25: Sydney 2024 Summer Temp - Classical vs. Quantum Quantum pdq (4, 1, 1) (3, 1, 1) (5, 1, 1) (6, 1, 1) (10, 1, 1) (7, 1, 1) (9, 1, 1) (8, 1, 1)

Table 26: Sydney 2024 Summer Classical VS Quantum Quantum pdq (3, 1, 1) (4, 1, 1) (5, 1, 1) (10, 1, 1) (6, 1, 1) (9, 1, 1) (8, 1, 1) (7, 1, 1)

DM Stats (MSE)

Sydney 2024 Summer Temp Classical VS Quantum DM Stats MSE blocks used loss DM stat p value classical mean loss quantum mean loss 160 MSE -3.41981 6.27E-04 8.1677 9.023377 161 MSE -2.943498 3.25E-03 8.177383 19.809067 161 MSE -2.456801 1.40E-02 8.177383 8.713587 161 MSE -1.544826 1.22E-01 8.177383 8.668534 155 MSE 0.610649 5.41E-01 8.230546 7.906464 161 MSE -0.577497 5.64E-01 8.177383 8.416235 159 MSE 0.555303 5.79E-01 8.258326 7.989777 161 MSE 0.175505 8.61E-01 8.177383 8.099677

DM Stats (MAE)

Sydney 2024 Summer Classical VS Quantum DM Stats MAE blocks used loss DM stat p value classical mean loss quantum mean loss 161 MAE -4.038151 5.40E-05 2.183719 2.875055 160 MAE -2.975823 2.92E-03 2.180798 2.280019 161 MAE -1.775833 7.58E-02 2.183719 2.245571 155 MAE 1.205682 2.28E-01 2.193426 2.106318 161 MAE -1.158140 2.47E-01 2.183719 2.237213 159 MAE 1.088950 2.76E-01 2.198211 2.124441 161 MAE 0.601969 5.47E-01 2.183719 2.145549 161 MAE -0.164042 8.70E-01 2.183719 2.193769

delta mean loss -0.855678 -11.631684 -0.536204 -0.491151 0.324082 -0.238851 0.268549 0.077707

delta mean loss -0.691336 -0.099221 -0.061852 0.087108 -0.053494 0.073770 0.038169 -0.010050

Record · ID 2610 · SHA-256 81187d337ac1746b
Conceptio Open Knowledge Archive — every document is proof-bundled with source, license, and retrieval metadata.