1 Introduction↩︎

Drought is among the most disruptive climate hazards, especially in semi arid regions such as South Australia, where it undermines agriculture, water security, and ecosystem resilience [1]. The Standardized Precipitation Evapotranspiration Index (SPEI) integrates precipitation with potential evapotranspiration and is widely used to characterise drought severity across time scales [2], [3].

Existing approaches face challenges in predicting SPEI one to three months in advance, especially in local regions where data is insufficient to train complex models. Traditional statistical models including ARIMA, SARIMA, and state space formulations [4] assume quasi stationarity and therefore struggle to capture the coexistence of short lived pulses and long term trends. Traditional machine learning algorithms such as support vector machines or random forests can model nonlinear relationships yet depend on hand crafted features and suffer in high dimensional, multivariate settings [5]. Deep architectures LSTMs [6], TCNs [7], and the growing family of Transformer variants [8], [9] can learn temporal dependencies end to end, but they require large datasets and are prone to overfitting in data sparse regions.

Large Time Series Foundation Models (TSFMs) have emerged as a promising alternative. Pre trained on millions of heterogeneous sequences, these models demonstrate strong zero-shot skill across finance, traffic, and meteorology [10]; recent examples include TimesFM [11] and Lag-Llama [12]. While physics aware encoder decoder designs such as GraphCast [13], FourCastNet [14], and Pangu-Weather [15] have advanced global numerical weather prediction, they require complex post processing pipelines to generate drought indices, making general purpose TSFMs more attractive for direct climate index forecasting.

Yet applying TSFMs to regional drought forecasting faces significant practical constraints that traditional adaptation approaches cannot address. Foundation model weights are often proprietary and inaccessible, preventing direct fine-tuning or architectural modifications. Regional climate records are typically short spanning decades rather than the centuries needed for robust model retraining while operational forecasting centers operate under strict computational budgets that preclude resource intensive adaptation methods. These constraints necessitate lightweight, inference time solutions that can enhance foundation model performance without requiring parameter access or extensive retraining.

Figure 1 illustrates these regional adaptation difficulties across three representative model families: a classical statistical model (SARIMA), a deep learning model (iTransformer), and a foundation model (TimesFM). The three differ because they rely on different inductive biases, yet all leave systematic rather than random residuals. The bias is most pronounced for SARIMA and iTransformer and smallest for the pre trained TimesFM, but even TimesFM retains a structured residual pattern. This indicates that regional temporal dynamics contain characteristics at different temporal resolutions that large scale pre training does not fully capture, motivating specialized inference time adaptation.

To address these practical constraints while preserving foundation models’ general capabilities, we introduce two lightweight inference time wrappers. The Stationarity aware multi-resolution Residual wrapper (SMR²) decomposes the regional forecasting task into multi-resolution temporal pattern learning, capturing stride specific regional dynamics that complement the global knowledge embedded in foundation models. MBB (Moving Block Bootstrap) constructs residual block perturbations of the target history and averages the corresponding frozen backbone forecasts.

In summary, our main contributions are as follows:

• We propose two plug-and-play inference time wrappers, SMR² for multi-resolution regional pattern learning and MBB for variance reducing block bootstrap ensembling, that require no access to the backbone foundation model weights.

• We present a theoretical analysis built on a single residual decomposition assumption: a convex risk bound with an explicit bias variance trade off that motivates the adaptive multi-resolution weighting, together with the standard block bootstrap variance identity that explains the ensemble’s stabilizing effect.

• We conduct comprehensive experiments on SPEI-30 forecasting across three South Australian sites, showing up to 26% mean squared error (MSE) reduction over the corresponding frozen backbone, consistently outperforming 14 statistical and deep learning baselines across three foundation model backbones.

2 Related Work↩︎

Drought Forecasting Approaches. Traditional statistical methods like ARIMA assume stationarity and struggle with nonlinear climate dynamics [1], while machine learning approaches (SVMs, Random Forests) require handcrafted features and lack temporal modeling capabilities [5]. Deep learning architectures including LSTMs [6], TCNs [7], and Transformers [8] enable end to end learning but require task specific retraining and suffer from limited regional data availability.

Time Series Foundation Models. Large scale pre trained models like TimesFM [11] and Lag-Llama [12] demonstrate strong zero-shot capabilities across domains. While domain specific weather models (GraphCast [13], FourCastNet [14]) excel at global numerical prediction, they generate gridded physical variables requiring extensive post processing for climate indices, introducing uncertainty sources and computational constraints that exceed regional system capabilities. This motivates our focus on general purpose TSFMs that can forecast climate indices directly, with regional adaptation handled outside the backbone.

Multi-resolution and Resampling Methods. multi-resolution analysis decomposes signals across temporal scales [16], with recent advances in multi scale attention [17] and residual frameworks like N-BEATS [18]. The Moving Block Bootstrap [19] preserves temporal dependencies when resampling serially correlated data [20]. These methods are therefore less suited to black-box backbone adaptation.

Our work builds on these ideas by applying multi-resolution residual correction and temporally aware resampling at inference time around a frozen backbone.

3 Preliminary↩︎

3.1 Problem Formulation↩︎

Let \(\{(\mathbf{c}_t, y_t)\}\) denote a multivariate climate time series, where each \(\mathbf{c}_t \in \mathbb{R}^{d}\) contains \(d=38\) predictor features (34 meteorological variables and 4 climate indices) and \(y_t \in \mathbb{R}\) is the scalar target SPEI value. For a time interval \(a{:}b\), \(y_{a:b}=(y_a,\ldots,y_b)\) denotes the corresponding target history. We define:

• \(L\): Length of the context window used for prediction

• \(L_{\text{train}}\): Reduced context length used during calibration phase (\(L_{\text{train}} < L\))

• \(\mathcal{T}_{\text{train}}, \mathcal{T}_{\text{val}}\): Index sets used for calibration and validation; held out test splits are used only for reporting

At each forecast time \(t\), we collect the context window \(\mathbf{X}_t = (\mathbf{c}_{t-L+1:t}, y_{t-L+1:t})\) and predict the next SPEI value with a frozen Time Series Foundation Model backbone \(f_{\theta}^{\text{base}}\), i.e., a fixed mapping \(\hat{y}_{t+1} = f_{\theta}^{\text{base}}(\mathbf{X}_t)\) from a context window to a forecast. We write \(\mathbf{X}_t^{(s)}\) for the same window after the stride \(s\) multi-resolution view construction is applied to every channel of \(\mathbf{X}_t\) (the \(d\) covariates and the target history), i.e.\(\mathbf{X}_t^{(s)} = (\mathbf{c}_{t-L+1:t}^{(s)}, y_{t-L+1:t}^{(s)})\), pairing the stride \(s\) covariate view \(\mathbf{c}_{t-L+1:t}^{(s)}\) with the stride \(s\) target history view \(y_{t-L+1:t}^{(s)}\), both constructed in the Method section below. We focus on one month ahead prediction.

4 Method↩︎

4.1 Framework Overview↩︎

To adapt frozen Time Series Foundation Models (TSFMs) for regional forecasting tasks, we propose two lightweight inference time wrappers that operate without modifying pre trained model parameters. Here, \(f_{\theta}^{\text{base}}\) denotes the frozen TSFM backbone (e.g., TimesFM, TimeGPT) that provides the base predictions; we write \(f_{\theta}^{\text{base}}\) throughout the Method to emphasize that it is not updated. The SMR² wrapper addresses regional pattern learning, while the MBB wrapper improves forecast robustness through block bootstrap ensembling. Figure 2 provides an overview of the two wrappers, highlighting their shared bagging style structure: diverse inference time series, a shared frozen backbone, and aggregation.

4.2 SMR² Wrapper↩︎

Foundation models trained on global datasets face adaptation difficulties when applied to regional forecasting tasks, owing to practical deployment constraints and the need to capture local temporal dynamics. Figure 3 demonstrates that regional time series exhibit non stationary shifts over time, different temporal strides reveal distinct distributional characteristics, and foundation model residuals show complementary error patterns across resolutions.

Instead of retraining the entire model, we can enhance regional forecasting by learning complementary temporal patterns across multiple resolutions that supplement the foundation model’s global knowledge. Different temporal resolutions capture different aspects of the underlying signal: fine grained patterns reveal short term fluctuations while coarse grained patterns expose long term trends. We decompose the regional enhancement problem into multi-resolution pattern learning, where each temporal stride contributes differently to the overall performance, and residual learning, where we learn a residual correction for each stride. We call the wrapper stationarity aware in this specific sense: the per segment mean and standard deviation used to build each stride view (Eqs. 2 –3 ) act as a local normalization that absorbs slow drifts in level and variance [21], rather than because it performs an explicit stationarity test, change point detection, or regime classification.

The SMR² approach addresses both aspects simultaneously. Let \(\mathcal{S}\) denote the set of temporal strides and let \(w^{(s)}\) be the ensemble weight assigned to stride \(s\): \[\label{eq:smr295forecast} \begin{align} \tilde{y}_{t+1}^{\mathrm{ens}} &= \sum_{s \in \mathcal{S}} w^{(s)}\, \tilde{y}_{t+1}^{(s)}, \\ \tilde{y}_{t+1}^{(s)} &= f_{\theta}^{\text{base}}(\mathbf{X}_t^{(s)}) + \hat{\epsilon}_{t+1}^{(s)}. \end{align}\tag{1}\] where \(\mathbf{X}_t^{(s)}\) is the stride \(s\) multi-resolution input view (of both the covariates and the target history), \(\tilde{y}_{t+1}^{(s)}\) is the corrected stride \(s\) forecast, \(\tilde{y}_{t+1}^{\mathrm{ens}}\) is the final SMR² ensemble forecast, and \(\hat{\epsilon}_{t+1}^{(s)}\) is the learned residual correction defined by the two phase procedure below.

4.2.0.1 multi-resolution Input Construction

We build the stride \(s\) input \(\mathbf{X}_t^{(s)}\) by applying multi-resolution downsampling to the full context window \(\mathbf{X}_t = (\mathbf{c}_{t-L+1:t}, y_{t-L+1:t})\), that is, to all \(d{+}1\) channels at once. Let \(\mathbf{x}_\tau = (\mathbf{c}_\tau, y_\tau) \in \mathbb{R}^{d+1}\) denote the input at a single time step \(\tau\), stacking the \(d\) covariates with the scalar target. For each stride \(s\), we partition the first \(s\lfloor L/s \rfloor\) positions of the window into \(\lfloor L/s \rfloor\) complete non overlapping segments. Any trailing remainder of length \(L-s\lfloor L/s \rfloor\) is not used to form an additional partial segment when computing the downsampled statistics. Writing \(I_j^{(s)} = \{(j-1)s+1, \ldots, js\}\) for the index set of segment \(j\), we compute, for each \(j \in \{1, 2, \ldots, \lfloor L/s \rfloor\}\), the per segment mean and standard deviation of the input: \[\begin{align} \bar{\mathbf{x}}_j^{(s)} &= \frac{1}{|I_j^{(s)}|} \sum_{k \in I_j^{(s)}} \mathbf{x}_{t-L+k} \tag{2} \\ \boldsymbol{\sigma}_j^{(s)} &= \sqrt{\frac{1}{|I_j^{(s)}|} \sum_{k \in I_j^{(s)}} \left(\mathbf{x}_{t-L+k} - \bar{\mathbf{x}}_j^{(s)}\right)^2} \tag{3} \end{align}\] Both operations act elementwise across the \(d{+}1\) channels, so the same stride \(s\) summary is applied to every covariate and to the target history alike.

Following downsampling, we perform temporal upsampling to obtain a length \(L\) representation. In the implementation, each segment’s statistics are repeated \(\lfloor L/\lfloor L/s\rfloor\rfloor\) times. If the repeated sequence is shorter than \(L\), the remaining positions are filled with the last complete segment statistic; if it is longer than \(L\), it is truncated. Thus, for the upsampled feature vector \(\mathbf{u}_{t-L+k}^{(s)}\), \[\begin{align} \mathbf{u}_{t-L+k}^{(s)} &= [\bar{\mathbf{x}}_{j}^{(s)},\boldsymbol{\sigma}_{j}^{(s)}], \\ j= & \min\left( \left\lceil \frac{k}{\lfloor L/\lfloor L/s\rfloor\rfloor} \right\rceil, \left\lfloor \frac{L}{s}\right\rfloor \right), \\ &\quad k=1,\ldots,L . \end{align} \label{eq:upsample}\tag{4}\] Concatenating mean with standard deviation provides richer geometric information than traditional delay embeddings, and stacking the upsampled vectors \(\mathbf{u}_{t-L+k}^{(s)}\) over \(k\) yields the multi-resolution input view \(\mathbf{X}_t^{(s)} \in \mathbb{R}^{L \times 2(d+1)}\) that is passed to the backbone. We write \(\mathbf{c}_{t-L+1:t}^{(s)} \in \mathbb{R}^{L \times 2d}\) and \(y_{t-L+1:t}^{(s)} \in \mathbb{R}^{L \times 2}\) for its covariate and target history sub blocks, so \(\mathbf{X}_t^{(s)} = (\mathbf{c}_{t-L+1:t}^{(s)}, y_{t-L+1:t}^{(s)})\).

4.2.0.2 Two Phase Prediction Framework

The two-phase framework separates residual learning from final forecasting. In the calibration phase, the frozen backbone is applied to each stride view to obtain a base forecast, and the difference between the observed target and this base forecast is stored as a stride-specific residual. In the inference phase, the same frozen backbone first produces the stride-specific base forecast \(\hat{y}^{(s)}_{t+1}\); a lightweight residual corrector then predicts \(\hat{\epsilon}^{(s)}_{t+1}\) from recent residuals and current stride-view summaries; the corrected stride forecast is finally obtained as \(\tilde{y}^{(s)}_{t+1}=\hat{y}^{(s)}_{t+1}+\hat{\epsilon}^{(s)}_{t+1}\). We use a reduced context length \(L_{\mathrm{train}}\) during calibration to compute initial residuals, while using the full context length \(L\) for final ensemble predictions.

During the calibration phase, we use the reduced context length \(L_{\text{train}}\) to compute calibration residuals for all \(t \in \mathcal{T}_{\text{train}}\). To avoid confusing this reduced window with the full inference window, let \(\mathbf{X}_{t-1,L_{\text{train}}}=(\mathbf{c}_{t-L_{\text{train}}:t-1}, y_{t-L_{\text{train}}:t-1})\) denote the length-\(L_{\text{train}}\) context ending at \(t-1\), and let \(\mathbf{X}_{t-1,L_{\text{train}}}^{(s)}\) be its stride \(s\) view: \[\epsilon_t^{(s)} = y_t - f_{\theta}^{\text{base}}(\mathbf{X}_{t-1,L_{\text{train}}}^{(s)})\] Here \(\epsilon_t^{(s)}\) is the realized base model residual at stride \(s\). Collecting these residuals gives the calibration residual set \(\mathcal{R}^{(s)}=\{\epsilon_t^{(s)}\}_{t\in\mathcal{T}_{\text{train}}}\). Calibration shares the forecast origin alignment of inference and uses no information from time \(t\) or later (no look ahead).

During the inference phase, for each stride \(s \in \mathcal{S}\), we run both models in parallel using the full context window \(L\): \[\begin{align} \hat{y}_{t+1}^{(s)} &= f_{\theta}^{\text{base}}(\mathbf{X}_t^{(s)}), \\ \hat{\epsilon}_{t+1}^{(s)} &= f^{\text{res}}\big(\boldsymbol{\phi}_t^{(s)}\big), \qquad \boldsymbol{\phi}_t^{(s)} = \big[\,\epsilon_{t-q+1:t}^{(s)},\; \psi(\mathbf{X}_{t}^{(s)})\,\big], \\ \tilde{y}_{t+1}^{(s)} &= \hat{y}_{t+1}^{(s)} + \hat{\epsilon}_{t+1}^{(s)}, \end{align}\] where the corrector input \(\boldsymbol{\phi}_t^{(s)}\) concatenates the \(q\) most recent stride \(s\) residuals \(\epsilon_{t-q+1:t}^{(s)}\) with summary statistics \(\psi(\cdot)\) of the stride \(s\) input view \(\mathbf{X}_t^{(s)}\) (the per channel means of its mean and standard deviation components, over both the covariates and the target history). The corrector \(f^{\text{res}}\) is a ridge regression with \(\ell_2\) penalty \(\lambda_{\text{ridge}}=0.1\), fit once on \(\mathcal{R}^{(s)}\) and then applied unchanged at inference, so the wrapper performs no test time gradient updates. The frozen backbone consumes the joint multi-resolution view \(\mathbf{X}_t^{(s)}\); for backbones whose interface accepts only a single stream (e.g., target history only or covariate only models), the unused component is simply omitted without otherwise changing the procedure.

4.2.0.3 Adaptive Ensemble Fusion

We choose weights that emphasize strides with low error. Our adaptive fusion strategy uses validation performance as a lightweight proxy for per stride reliability. Let \(R^2(s)\) and \(\mathrm{MSE}^{(s)}\) denote the validation coefficient of determination and mean squared error for stride \(s\), and let \(\varepsilon>0\) be a small stabilizer. The choice of \(R^2\)-based weighting reflects that each stride’s contribution should be proportional to its explained variance, naturally emphasizing strides that better capture the underlying signal structure. When systematic biases dominate across all strides, every \(R^2(s)\) becomes non positive (each stride does worse than predicting the mean), so the \(R^2\)-based weights are undefined or negative. We therefore fall back to inverse MSE weighting, which is justified because it remains strictly positive and normalisable, still assigns larger weight to lower error strides, and recovers a sensible ranking precisely in the regime where \(R^2\) is uninformative: \[\label{eq:r2weights} w^{(s)} = \begin{cases} \dfrac{\max(R^2(s),0)}{\sum_{s' \in \mathcal{S}} \max(R^2(s'),0)} & \text{if } \exists\, s: R^2(s) > 0 \\[8pt] \dfrac{1/(\mathrm{MSE}^{(s)} + \varepsilon)}{\sum_{s' \in \mathcal{S}} 1/(\mathrm{MSE}^{(s')} + \varepsilon)} & \text{otherwise} \end{cases}\tag{5}\] The \(\max(\cdot,0)\) truncation guarantees that every weight is non negative and that the weights form a valid simplex, so no stride can contribute a negative coefficient. The final ensemble output then uses these adaptive weights in Eq. 1 .

4.3 MBB Wrapper↩︎

Foundation model forecasts can be sensitive to the exact target history supplied at inference. MBB reduces this sensitivity by constructing \(K\) temporally coherent perturbations of the target history, feeding each perturbed history into the same frozen backbone, and averaging the resulting forecasts. The covariate history is kept unchanged; only the target history is perturbed.

Figure 4 illustrates this perturb and average workflow: block-resampled residual sequences generate multiple perturbed target histories, each perturbed history induces a frozen-backbone forecast, and their average forms the final MBB forecast.

The perturbations are generated from validation residuals of the frozen backbone. Specifically, we first collect residuals \(\{r_t = y_t-\hat{y}_t\}_{t\in\mathcal{T}_{\mathrm{val}}}\), where \(\hat{y}_t\) is the corresponding one-step forecast from the frozen backbone. Since these residuals are serially correlated, sampling individual residuals independently would break local autocorrelation. MBB therefore samples short consecutive residual blocks rather than individual residuals. Sampling and concatenating such blocks produces residual perturbation sequences that preserve local temporal dependence while generating diverse target-history variants.

4.3.0.1 Residual Block Construction

Following standard MBB practice [19], [22], we set the block length to \(b = \lfloor n^{1/3} \rfloor\) where \(n = |\mathcal{T}_{\text{val}}|\), a common asymptotic rule of thumb under standard dependence assumptions. We construct overlapping candidate blocks: \[B_i = \{r_i, r_{i+1}, \ldots, r_{i+b-1}\}, \quad i \in \{1, \ldots, n-b+1\}\] This overlapping construction lets the finite validation residual series contribute many possible local error patterns.

4.3.0.2 Bootstrap Resampling Strategy

To generate a perturbed sequence of length \(L\), we sample \(\lceil L/b \rceil\) blocks with replacement. The sampling process creates the ordered block sequence \[\mathbf{B}^{(k)} = (B_1^{(k)}, B_2^{(k)}, \ldots, B_{\lceil L/b \rceil}^{(k)})\] where each \(B_j^{(k)}\) is randomly selected from the available blocks. These blocks are then concatenated and truncated to form residual sequences of the required length: \[r_{1:L}^{(k)} = [B_1^{(k)}, B_2^{(k)}, \ldots]_{1:L}\] where \([\,\cdot\,]\) denotes the same vector concatenation used earlier (e.g., in the multi-resolution input construction) and \([\,\cdot\,]_{1:L}\) truncates the concatenated blocks to length \(L\). This process generates \(K\) different residual sequences used to form the perturbed histories below.

4.3.0.3 Ensemble Prediction

For each bootstrap sample \(k\), we construct perturbed target histories by adding the position indexed resampled residual sequence to the original predictions over the history window: \[\tilde{y}_{t-L+i}^{(k)} = \hat{y}_{t-L+i} + r_i^{(k)}, \quad i = 1, \ldots, L\] Stacking these over the window gives the perturbed target history \(\tilde{y}_{t-L+1:t}^{(k)} = (\tilde{y}_{t-L+1}^{(k)}, \ldots, \tilde{y}_{t}^{(k)})\), that is, the foundation model’s own predicted history corrected by the \(k\)th resampled residual sequence. Here \(r_i^{(k)}\) is the \(i\)th element of \(r_{1:L}^{(k)}\), and \(\hat{y}_{t-L+i}\) are the frozen backbone’s rolling one step predictions over the history window, obtained under the same one step protocol as inference. Pairing it with the covariates gives the \(k\)th perturbed input \(\tilde{\mathbf{X}}_t^{(k)} = (\mathbf{c}_{t-L+1:t}, \tilde{y}_{t-L+1:t}^{(k)})\), which is fed to the frozen foundation model: \[\begin{align} \hat{y}_{t+1}^{(k)} &= f_{\theta}^{\text{base}}(\tilde{\mathbf{X}}_t^{(k)}) \\ \hat{y}_{t+1}^{\mathrm{MBB}} &= \frac{1}{K} \sum_{k=1}^{K} \hat{y}_{t+1}^{(k)} \end{align}\] The second line is the ensemble point forecast; the variance argument for this average is stated in the theoretical analysis below.

4.4 Theoretical Analysis↩︎

Because the wrappers are deliberately lightweight, the analysis below is to make explicit why a black-box correction of a frozen TSFM should help: why the residual of a frozen backbone is learnable, why combining multiple temporal resolutions is sound, and why a block bootstrap ensemble stabilizes the point forecast. Everything rests on a single assumption and two short, assumption light results; the derivations are in Appendix A.

4.4.0.1 Residual notation in the analysis.

We use different residual symbols because the two wrappers use residuals for different purposes. In SMR\(^2\), \(\epsilon_t^{(s)}\) denotes the reduced-context, stride-specific calibration residual defined in Eq. (5). It is used to fit the stride-specific residual corrector and to produce the learned correction \(\hat{\epsilon}_{t+1}^{(s)}\) in Eq. (7). In MBB, \(r_t\) denotes the full-context validation residual of the frozen backbone, used only for constructing temporally coherent block-bootstrap perturbations. Thus, \(\epsilon_t^{(s)}\) is a supervised correction target, whereas \(r_t\) is a resampling object.

This is consistent with what Figures 1 and 3 display: every backbone in Figure 1 (SARIMA, iTransformer, TimesFM) leaves a systematic, non random residual, and Figure 3 shows the residual distribution changing shape across strides, so the multi-resolution components carry complementary information rather than copies of the same signal.

4.4.0.2 SMR²: a convex risk bound.

Recall the per stride corrected forecast \(\tilde{y}_{t+1}^{(s)} = \hat{y}_{t+1}^{(s)} + \hat{\epsilon}_{t+1}^{(s)}\) and the ensemble \(\tilde{y}_{t+1}^{\mathrm{ens}} = \sum_{s\in\mathcal{S}} w^{(s)}\, \tilde{y}_{t+1}^{(s)}\) from the SMR² wrapper above, with non negative weights obeying \(\sum_{s\in\mathcal{S}} w^{(s)} = 1\).

Inequality ?? (a direct consequence of Jensen’s inequality, since the weights lie on the simplex) guarantees the ensemble risk is controlled by the weighted average of the per stride risks. The decomposition ?? motivates allocating mass to low error strides while accounting for inter stride dependence: when the covariance structure is favourable, mixing weakly correlated reliable strides can reduce the variance term \(\mathbf{w}^{\!\top}\boldsymbol{\Sigma}\,\mathbf{w}\) relative to a single chosen stride. The \(R^2\)/inverse MSE weights in Eq. 5 serve as a lightweight validation proxy for this allocation by favouring low risk resolutions, while remaining covariance agnostic. Intuitively, coarser strides average out high frequency fluctuations and expose slower regional trends, so different resolutions correct different parts of the residual in Eq. ?? .

4.4.0.3 MBB: variance stabilization.

The MBB wrapper averages \(K\) forecasts obtained from block bootstrapped residual perturbations of the target history. For any such ensemble, \[\label{eq:mbbvar} \mathrm{Var}\!\left(\tfrac{1}{K}\sum_{k=1}^{K}\hat{y}_{t+1}^{(k)}\right) = \frac{1}{K^2}\sum_{i,j}\mathrm{Cov}\big(\hat{y}_{t+1}^{(i)}, \hat{y}_{t+1}^{(j)}\big),\tag{6}\] which equals \(\tfrac{1}{K}\mathrm{Var}(\hat{y}_{t+1})\) when the bootstrap forecasts are independent and, for approximately exchangeable forecasts of comparable variance, falls below the single forecast variance whenever their pairwise correlation is below one. Under standard stationarity and mixing conditions the Moving Block Bootstrap is a consistent resampling scheme for dependent series [19], [22], [23]; preserving within block autocorrelation (block length \(b=\lfloor n^{1/3}\rfloor\)) keeps the perturbed histories realistic, so averaging over them stabilizes the point forecast without touching the backbone.

5 Experiments↩︎

5.1 Study Area↩︎

We focus our evaluation on South Australia, shown in Figure 6, a region characterized by high drought vulnerability and complex climate variability driven by multiple interacting modes (ENSO, IOD). South Australia occupies a unique position in Australia’s climate system, experiencing a temperate climate with dry summers in the south transitioning to a persistently dry hot desert climate in the north [24]. The region’s climate exhibits high spatial and temporal variability, with documented warming trends since the 1970s and strong dependence on multiple large scale climate drivers. It provides high quality long term meteorological data and represents a critical water security area where improved drought forecasting has significant practical value, while recent intensifying drought conditions create challenging adaptation scenarios for globally pre trained foundation models. All datasets used in this study are publicly available and reproducible through standardized preprocessing pipelines, enabling transparent evaluation and potential deployment across other semi arid regions.

5.2 Experimental Setup↩︎

Dataset. We evaluate on monthly SPEI-30 prediction across 3 South Australian locations over 1982–2018. The target is the SPEI-30 variable from the global daily SPEI-GD database [25] (0.25\(^\circ\) spatial resolution, natively spanning 1982–2021), converted to a monthly cadence and restricted to 1982–2018. The dataset comprises 39 variables: the target SPEI plus \(d=38\) predictors, namely 34 meteorological variables from the NCEP–NCAR Reanalysis 1 [26] and 4 climate indices (e.g., ENSO and IOD [27], [28]). We use a chronological 70/10/20 train/validation/test split; each forecast uses a context window of length \(L\) covering the available training history.

Baselines. We compare our wrappers against three categories of models: (i) classical statistical and recurrent baselines (ARIMA [4], ETS, LSTM [6], GRU, and TCN [7]); (ii) Transformer based architectures including standard Transformer [8], Autoformer [29], Crossformer [30], DLinear [31], FiLM [32], iTransformer [33], PatchTST [34], TimesNet [9], and TSMixer [35]; and (iii) foundation models including TabPFN [36], TimeGPT [37], and TimesFM [11]. All models are evaluated using identical data splits and preprocessing pipelines.

Implementation Details. SMR² uses strides \(\mathcal{S} = \{1, 2, 4, 6, 8\}\) selected to capture multi-resolution temporal patterns from monthly to seasonal dynamics (see Appendix B for details). The ridge residual corrector \(f^{\text{res}}\) uses penalty \(\lambda_{\text{ridge}} = 0.1\), residual lag length \(q = 6\), and stabilizer \(\varepsilon = 10^{-6}\). MBB sets block length \(b = \lfloor n^{1/3} \rfloor\) and \(K = 100\) bootstrap samples. Both wrappers add only inference time computation; runtime, memory, hardware, and software details are reported in Appendices E–F.

Evaluation Metrics. We report Mean Squared Error (MSE), Mean Absolute Error (MAE), and coefficient of determination (R\(^2\)); the evaluation protocol is detailed in Appendix D.

5.3 Main Experiments↩︎

Forecasting Strategy. We use a rolling one-step-ahead forecasting approach where predictions are made sequentially across the test period, with the context window sliding forward after each prediction to evaluate performance throughout the entire test set.

For baseline models, we train each architecture from scratch on the training set using their original best configurations. Traditional deep learning models (LSTM, GRU, TCN) are trained for 100 epochs with early stopping based on validation performance. Transformer based models use the same training protocol with architecture specific hyperparameters as reported in their original papers. Foundation models (TabPFN-ts V1.0, TimeGPT-1, and TimesFM v1.0) are used in their pre trained form without additional fine-tuning. This is not meant to understate their capacity: fine-tuning is precisely what the deployment setting we target rules out, since the backbone weights are often proprietary, the regional records are too short for stable fine-tuning, and operational compute budgets are limited. Evaluating the foundation models zero-shot therefore reflects the realistic regime in which our wrappers would be used, and provides a controlled comparison in which both the baselines and our method leave the backbone frozen; the wrappers add only inference time computation on top of the same frozen predictions. For each foundation model, we evaluate three configurations: baseline without modification, with SMR² wrapper, and with MBB wrapper applied at inference time. We evaluate the two wrappers separately to isolate the effect of each; they are alternative lightweight adaptations rather than a combined pipeline.

We evaluate performance across three representative locations in South Australia that capture different climate zones. For each location, we compute all metrics separately and report location specific results in Table 1. Figure 5 visualizes the predicted versus observed SPEI at Location 1 and provides representative visual evidence consistent with the quantitative improvements in Table 1.

Table 1 shows that frozen foundation models already outperform most task trained baselines, while SMR² and MBB further improve every frozen backbone across all three sites, indicating complementary regional corrections without parameter updates.

5.4 Ablation Study↩︎

Table 2 presents an ablation study on TimesFM v1.0 at Location 1, systematically evaluating each component’s contribution. For SMR², incremental addition of components shows progressive improvements: single stride residual correction (S-1) reduces MSE by 2.0%, S-2 uses multiple strides with uniform weights and achieves 4.1% reduction, and the full configuration (S-Full), which combines adaptive \(R^2\) weighting, two phase calibration, and the mean/std multi-resolution embedding, delivers 19.2% MSE reduction. For MBB, the ablation isolates the effect of block sampling: IID bootstrap (B-1) provides minimal improvement (1.8%), while the full MBB configuration (B-Full) with overlapping resampling and the standard block length \(b=\lfloor n^{1/3}\rfloor\) achieves 19.9% MSE reduction. Both wrappers show substantial and complementary improvements over the frozen baseline.

5.5 Generalization to Other Regions↩︎

To evaluate cross regional robustness, we further test our inference time wrappers with the TimesFM v1.0 backbone on three climatically distinct regions: Indonesia (tropical monsoon), North Africa (arid desert), and the U.S. West Coast (Mediterranean). These regions exhibit markedly different rainfall seasonality and drought dynamics, providing a stringent generalization test.

All regions follow the same experimental protocol as South Australia, including identical predictor target construction, preprocessing, rolling origin evaluation, and hyperparameter settings. Monthly reanalysis variables are used as predictors, with the SPEI-30 variable from the global daily SPEI-GD dataset (0.25°, converted to a monthly cadence over 1982–2018) as the target. The backbone TimesFM model remains fully frozen.

As summarized in Table 3, SMR² reduces MSE in every region, by roughly 16 to 41%, while MBB also improves the frozen backbone under the same cross regional protocol. The largest relative gains appear on the U.S. West Coast and in North Africa, with a smaller but still clear improvement on the tropical monsoon Indonesia case. These results demonstrate that the proposed wrappers generalize effectively across heterogeneous climate regimes without fine-tuning or parameter updates.

6 Social Impact↩︎

This work supports the responsible use of pre trained time series foundation models for regional drought analysis under realistic data and deployment constraints.

Improved forecast reliability facilitates more informed interpretation of model outputs in downstream applications such as water resource management, agricultural planning, and disaster preparedness, helping practitioners assess risk under climate variability.The framework may also help extend scientific foundation models to regional settings where intrusive model modification is infeasible.

7 Limitations, Ethical Considerations, and Conclusion↩︎

This work evaluates monthly SPEI forecasting across multiple regions and demonstrates stable cross regional generalization within the evaluated scope. The study focuses on a limited set of temporal resolutions and climate settings, and performance under substantially different temporal granularities, climate regimes, or hydrological targets has not been systematically validated. In addition, the framework relies on sufficient historical observations for inference time calibration, and may be affected by data scarcity or by inherent biases of frozen foundation models under extreme regime shifts.

All data used are publicly available and do not involve human subjects or personal information. As drought forecasts may inform high stakes decisions such as water resource management and agricultural planning, predictions should be used as decision support in conjunction with local observations and domain expertise.

In summary, we proposed a lightweight inference time adaptation framework for applying frozen time series foundation models to regional drought forecasting. Without modifying model parameters or training objectives, the SMR² and MBB wrappers effectively improve predictive accuracy in data scarce regional climate settings. Owing to its modular, backbone agnostic design, the framework is broadly applicable to climate science tasks that rely on frozen or black-box models.

GenAI Disclosure↩︎

We used a large language model (LLM) solely as a writing assistant to translate and improve grammar, fluency, and consistency of the manuscript text. All ideas, methods, analyses, and conclusions are by the authors; all citations and numerical values were manually verified against primary sources. No non-public data were provided to the tool.

Acknowledgments↩︎

This work was supported by the ARC Discovery Project DP230101122 and the Adelaide University Research Training Program (RTP) Scholarship. We gratefully acknowledge the continued support from the CSIRO Environment and Technology Research Units.

Appendix A: Residual Decomposition and Risk Bound↩︎

This appendix states the assumption and the two short results used in the main text, and proves the only non trivial one (the SMR² risk bound). The Moving Block Bootstrap properties are classical and are cited rather than reproved.

Assumption 1 (multi-resolution residual decomposition). The one step backbone residual \(r_{t+1}=y_{t+1}-f_{\theta}^{\text{base}}(\mathbf{X}_t)\) decomposes as \(r_{t+1}=\sum_{s\in\mathcal{S}} r_{t+1}^{(s)}+\eta_{t+1}\), with stride specific predictable components \(r_{t+1}^{(s)}\) and irreducible noise \(\eta_{t+1}\).

Proposition 1 (Ensemble risk bound). Let \(e^{(s)}=\tilde{y}_{t+1}^{(s)}-y_{t+1}\) and let \(\mathbf{w}=(w^{(s)})_{s\in\mathcal{S}}\) satisfy \(w^{(s)}\ge0\), \(\sum_s w^{(s)}=1\). Then \[\begin{gather} \mathbb{E}\big[(\tilde{y}_{t+1}^{\mathrm{ens}}-y_{t+1})^2\big]\le \sum_{s\in\mathcal{S}} w^{(s)}\,\mathbb{E}\big[(e^{(s)})^2\big],\\[2pt] \mathcal{L}_{\mathrm{ens}}=(\mathbf{w}^{\!\top}\boldsymbol{\beta})^2+\mathbf{w}^{\!\top}\boldsymbol{\Sigma}\,\mathbf{w}. \end{gather}\] where \(\boldsymbol{\beta}=(\mathbb{E}[e^{(s)}])_{s}\) and \(\boldsymbol{\Sigma}=\mathrm{Cov}\big((e^{(s)})_{s}\big)\).

Proof. Because the weights lie on the simplex, the ensemble error \(\tilde{y}_{t+1}^{\mathrm{ens}}-y_{t+1}=\sum_{s} w^{(s)} e^{(s)}=\mathbf{w}^{\!\top}\mathbf{e}\) is a convex combination of the per stride errors \(\mathbf{e}=(e^{(s)})_{s}\). Since \(x\mapsto x^2\) is convex, Jensen’s inequality gives \(\big(\sum_{s} w^{(s)} e^{(s)}\big)^2\le \sum_{s} w^{(s)} (e^{(s)})^2\) pointwise; taking expectations yields the first inequality. For the identity, \[\mathcal{L}_{\mathrm{ens}}=\mathbb{E}\big[(\mathbf{w}^{\!\top}\mathbf{e})^2\big]=\big(\mathbb{E}[\mathbf{w}^{\!\top}\mathbf{e}]\big)^2+\mathrm{Var}(\mathbf{w}^{\!\top}\mathbf{e})=(\mathbf{w}^{\!\top}\boldsymbol{\beta})^2+\mathbf{w}^{\!\top}\boldsymbol{\Sigma}\,\mathbf{w}.\] When the covariance structure is favourable, for example when reliable low variance strides have sufficiently weak positive or negative correlation, suitable simplex weights can satisfy \(\mathbf{w}^{\!\top}\boldsymbol{\Sigma}\,\mathbf{w}<\min_{s\in\mathcal{A}}\Sigma_{ss}\) for a subset \(\mathcal{A}\) of reliable strides. In other cases, the Jensen bound above still guarantees that the ensemble risk is controlled by the weighted average of the per stride risks. \(\square\)

MBB variance stabilization (standard result). The MBB forecast is the average \(\hat{y}_{t+1}^{\mathrm{MBB}}=\tfrac{1}{K}\sum_{k=1}^{K}\hat{y}_{t+1}^{(k)}\) of forecasts from block bootstrapped residual perturbations, with variance given by the elementary identity \[\mathrm{Var}\big(\hat{y}_{t+1}^{\mathrm{MBB}}\big)=\frac{1}{K^2}\sum_{i,j}\mathrm{Cov}\big(\hat{y}_{t+1}^{(i)},\hat{y}_{t+1}^{(j)}\big),\] which equals \(\tfrac{1}{K}\mathrm{Var}(\hat{y}_{t+1})\) for independent replicates and, for approximately exchangeable replicates of comparable variance, stays below the single forecast variance whenever the pairwise correlation is below one. Consistency of the Moving Block Bootstrap for stationary, mixing sequences is a classical result of Künsch [19], Lahiri [22], and Politis and Romano [23]; block length calculations often yield order \(b=O(n^{1/3})\) under specific objectives and dependence assumptions. We use the scheme as a black-box, input perturbation ensemble for a frozen TSFM rather than reproving its properties.

Appendix B: Algorithms and Implementation Details↩︎

Algorithms 7 and 8 use the same stride view construction as Eqs. 2 –4 : for the current window length \(m\), the construction pools only the first \(s\lfloor m/s\rfloor\) positions into complete stride \(s\) blocks over all \(d{+}1\) channels, repeats each block statistic \(\lfloor m/\lfloor m/s\rfloor\rfloor\) times, and forward fills or truncates to return a length \(m\) view.

Algorithm 1 SMR² Residual Extraction (Concise)

Algorithm 2 SMR² Inference

Algorithm 3 MBB Wrapper

Appendix C: Hyperparameter Configuration↩︎

We summarize below the key hyperparameters used for SMR² and MBB. SMR² uses strides \(\mathcal{S}=\{1,2,4,6,8\}\) aligned with dominant spectral components of monthly SPEI. Residual correction applies ridge regression (\(\lambda_{\text{ridge}}=0.1\), residual lag length \(q=6\)) on a feature vector of recent residuals and stride \(s\) input view summary statistics, with normalization \(\varepsilon=10^{-6}\), and ensemble weights based on validation \(R^2(s)\). MBB adopts block length \(b=\lfloor n^{1/3} \rfloor\) and \(K=100\) bootstrap samples, averaged into the ensemble point forecast. All models are frozen during inference, and hyperparameters are tuned only on validation data.

Appendix D: Dataset and Evaluation Protocol↩︎

Data sources. All datasets used in this study are publicly available. Meteorological predictors are obtained from the NCEP–NCAR Reanalysis 1 dataset. Drought targets use the SPEI-30 variable from the global daily SPEI-GD database (0.25\(^\circ\), native span 1982–2021), converted to monthly cadence. The study region covers South Australia, with representative grid points at \((-26.125^\circ,129.125^\circ)\), \((-29.125^\circ,134.875^\circ)\), and \((-35.625^\circ,138.875^\circ)\). The SPEI product is accessible at https://doi.org/10.5281/zenodo.8060268. Documentation and downloads for NCEP–NCAR Reanalysis 1 are available from NOAA at https://www.psl.noaa.gov/data/gridded/data.ncep.reanalysis.html. Public climate indices (e.g., ENSO, IOD) are provided by operational centers such as NOAA CPC.

Temporal coverage. Monthly data from 1982–2018, using the SPEI-30 variable from the daily SPEI-GD product, which natively spans 1982–2021.

Preprocessing. Missing values are linearly interpolated and predictors standardized. The target is the SPEI-30 series from SPEI-GD, converted to the monthly cadence used for evaluation; no additional target smoothing is applied beyond this cadence conversion.

Evaluation protocol. All models are evaluated with rolling one-step-ahead forecasting over the held out test period. For all Deep Learning methods implementation, the THUML Time-Series-Library [38] provides consistent data loaders, training scripts, and reference configurations that cover most baselines. For our method implementation, we provide it in our code repo.

Appendix E: Computational Cost and Deployment Notes↩︎

Runtime and hardware. Experiments are run on a single workstation equipped with an NVIDIA GeForce RTX 3090 GPU with 24 GB of memory, an Intel 12700K CPU, 64 GB of DDR4 3200 memory, and a 2 TB NVMe solid state drive. The operating system is Ubuntu 20.04.3 LTS. Average inference latency is below 3 s per site, with memory usage under 4 GB.

Energy and sustainability. Because SMR² and MBB operate exclusively at inference time, they avoid the repeated forward and backward passes over many epochs that fine-tuning requires, and instead add only a small number of additional forward passes (the \(K\) bootstrap evaluations for MBB and the \(|\mathcal{S}|\) stride evaluations for SMR²) to a single frozen model inference. Estimating the cost of fine-tuning as \(2\times\) the forward cost per step (the backward pass) over the typical number of update steps used for the deep learning baselines shows that the wrappers avoid the dominant training time computation; the exact energy saving depends on the fine-tuning budget and hardware. We therefore report measured latency and memory, and treat energy reduction as a qualitative deployment benefit rather than a measured percentage.

Deployment. Both wrappers are implemented as lightweight, plug-and-play modules for the frozen foundation model backbones evaluated in this paper (TabPFN, TimeGPT, and TimesFM), and can run efficiently on CPU only systems.

Appendix F: Reproducibility Checklist↩︎

Code and data availability. All scripts, processed data, and documentation are available at:

https://github.com/Wentao-Gao/Lightweight_Wrappers

Software environment. The software stack uses Python 3.10, PyTorch 2.0.1 with CUDA 11.8, NumPy 1.24.3, Pandas 2.0.2, and Scikit learn 1.2.2. Random seeds are fixed for NumPy and PyTorch so that runs are deterministic where supported. The seed affects trainable baselines and MBB resampling; frozen foundation model backbones are deterministic at inference.

Reproducibility scripts. The code repository includes configuration files and single command scripts that reproduce all scores and figures reported in the paper without retuning.

References↩︎