PAR(p) Inflow Model

Purpose

This spec defines the Periodic Autoregressive model of order $p$ (PAR(p)) used to capture temporal correlation in inflow time series, including the model definition, parameter semantics, the relationship between stored and computed quantities, the fitting procedure, model order selection, and validation invariants. Section 9 describes the optional PAR(p)-A extension that adds a single annual coefficient on top of the periodic AR structure to capture multi-year hydrological persistence.

1. Model Definition

The Periodic Autoregressive model of order p (PAR(p)) captures temporal correlation in inflow time series while accounting for seasonal variation in parameters. For hydro $h$ at stage $t$ corresponding to season $m(t)$:

$$a_{h,t} = \mu_{m(t)} + \sum_{\ell=1}^{p} \psi_{m(t),\ell} \left( a_{h,t-\ell} - \mu_{m(t-\ell)} \right) + \sigma_{m(t)} \cdot \varepsilon_t$$

where:

• $a_{h,t}$: Incremental inflow at stage $t$ (m³/s)
• $\mu_{m(t)}$: Seasonal mean for season $m(t)$
• $\psi_{m(t),\ell}$: Autoregressive coefficient for lag $\ell$ in season $m(t)$
• $\sigma_{m(t)}$: Innovation standard deviation for season $m(t)$ (derived at load — see section 3)
• $\varepsilon_t \sim \mathcal{N}(0, 1)$: Innovation (standardized noise)
• $m(t)$: Season index for stage $t$ (e.g., month 1–12)

The model order $p$ can vary by season and by hydro plant.

2. Parameter Set

For each hydro $h$ and each season $m \in \{1, \ldots, M\}$ (e.g., $M = 12$ for monthly, $M = 52$ for weekly), the complete PAR(p) model requires:

Parameter	Symbol	Description
Seasonal mean	$\mu_m$	Mean inflow for season $m$
AR coefficients	$\psi_{m,1}, \ldots, \psi_{m,p}$	Autoregressive coefficients
Innovation standard deviation	$\sigma_m$	Scale of the innovation term — derived, not independent (section 3)

3. Stored vs. Computed Quantities

The data model stores seasonal sample statistics and standardized AR coefficients — nothing else. The innovation scale is not a third stored quantity: the coefficient file holds the scale-invariant $\psi^*$ alone; at load the dimensionless innovation scale $r_m$ is derived from $\psi^*$ by a periodic-ACF closure (see “Deriving the innovation scale” below), then everything is converted to original-unit $\psi$ and $\sigma$ using the seasonal stats and consumed by the LP stage subproblem.

Two planes. Every PAR($p$) parameter sits on one of two planes, and the split is deliberate:

• Conditioning plane — the seasonal mean $\mu_m$ and seasonal sample standard deviation $s_m$ (both m³/s): the level and magnitude of the series. A study may re-condition these — e.g. a climate scenario that shifts both the seasonal mean and the seasonal variability.
• Dynamics plane — the standardized coefficients $\psi^*_{m,\ell}$ and the standardized innovation scale $r_m = \sigma_m/s_m$ (both dimensionless): the shape of the temporal dependence, independent of magnitude. Only $\psi^*_{m,\ell}$ is an independent parameter here — $r_m$ is pinned by $\psi^*_{m,\ell}$ itself.

Runtime re-couples the two via $\psi_{m,\ell} = \psi^*_{m,\ell}\,s_m/s_{m-\ell}$ and $\sigma_m = s_m\,r_m$, so re-conditioning the first plane rescales both the coefficients and the noise while preserving the correlation structure. The standardization basis is the seasonal std: $\psi^*$ is the coefficient of the process normalised by $s_m$ — the classical standardization of the periodic AR literature, and the same basis an externally-fitted model supplies (see the Inputs & Outputs tab).

Stored in input files

These are provided in inflow_seasonal_stats.parquet and inflow_ar_coefficients.parquet:

Stored quantity	Column	File	Symbol	Description
Seasonal sample mean	mean_m3s	inflow_seasonal_stats	$\mu_m = \bar{a}_m$	Mean of historical observations for season $m$
Seasonal sample std	std_m3s	inflow_seasonal_stats	$s_m$	Standard deviation of historical observations for season $m$
AR coefficients	coefficient	inflow_ar_coefficients	$\psi^*_{m,\ell}$	AR coefficient standardized by seasonal std — the direct Yule-Walker output

The AR order $p_m$ is not stored explicitly. It is derived at runtime from the count of coefficient rows per (hydro_id, stage_id) group in inflow_ar_coefficients.parquet.

The standardized coefficient $\psi^*_{m,\ell}$ is the direct output of the Yule-Walker fitting procedure (see section 5.4). It is dimensionless — the coefficient of the standardized process $(a_{h,t} - \mu_m) / s_m$. The relationship to the original-unit coefficient $\psi_{m,\ell}$ used in the LP is:

$$\psi_{m,\ell} = \psi^*_{m,\ell} \cdot \frac{s_m}{s_{m-\ell}}$$

Deriving the innovation scale

Unlike $\mu_m$, $s_m$, and $\psi^*_{m,\ell}$, the innovation scale $r_m$ is not read from a file. It is pinned by $\psi^*$ itself under the model’s unit-marginal-variance contract — every season’s standardized process $(a_{h,t}-\mu_m)/s_m$ has unit variance — via a periodic-ACF closure.

Extend the recursion that defines $\psi^*_{m,\ell}$ (section 5.4) to every lag $\ell \ge 1$, not just $\ell \le p_m$:

$$\rho_m(\ell) = \sum_{j=1}^{p_m} \psi^*_{m,j}\, \rho_{(m-\min(j,\ell)) \bmod M}\bigl(|\ell - j|\bigr), \qquad \rho_{m'}(0) = 1 \text{ for every season } m'$$

Solving this system jointly across every season $m = 0,\ldots,M-1$ and every lag $\ell = 1,\ldots,p_{\max}$ (where $p_{\max} = \max_m p_m$) yields the model’s own implied periodic autocorrelation function $\rho_m(\ell)$ — the autocorrelations the process would have if $\psi^*$ described an exactly stationary periodic AR process. The innovation scale follows from the same variance decomposition used in the fitting step (section 5.5):

$$r_m^2 = 1 - \sum_{\ell=1}^{p_m} \psi^*_{m,\ell}\, \rho_m(\ell)$$

A season with $p_m = 0$ contributes no terms to the sum, so $r_m = 1$ — a white-noise season gets no variance reduction from an AR part it does not have.

Why the innovation scale is derived, not stored

Under the unit-marginal-variance contract, $r_m$ is not a free parameter of the model — it is exactly the value that makes the seasonal marginal variance implied by $\psi^*$ equal $s_m^2$. Storing $r_m$ independently of $\psi^*$ would let the two drift out of consistency (e.g. a hand-supplied ratio that does not match the accompanying coefficients); deriving it instead guarantees the process’s implied marginal std in every season equals the stored $s_m$ exactly, whatever $\psi^*$ is.

Computed at runtime

From the stored quantities and the derived $r_m$, the LP requires two additional quantities computed once at initialization:

Original-unit AR coefficients (for LP constraint matrix entries):

$$\psi_{m,\ell} = \psi^*_{m,\ell} \cdot \frac{s_m}{s_{m-\ell}}$$

Innovation standard deviation (the noise scale):

$$\sigma_m = s_m \cdot r_m$$

No further autocorrelation values are needed beyond the closure derivation above. All required quantities — $\psi_{m,\ell}$, $\sigma_m$, and the $r_m$ they depend on — are derived solely from the stored seasonal stats and AR coefficient file, with no historical data and no separate noise-scale input.

LP coefficients

The stored standardized coefficients $\psi^*_{m,\ell}$ are converted to original-unit $\psi_{m,\ell}$ at runtime (see section 7.2), and these enter the LP directly (see LP Formulation). The LP equation is:

$$a_h = \underbrace{\left( \mu_m - \sum_{\ell=1}^{p} \psi_{m,\ell} \mu_{m-\ell} \right)}_{\text{deterministic base}} + \underbrace{\sum_{\ell=1}^{p} \psi_{m,\ell} \cdot a_{h,\ell}}_{\text{lag contribution}} + \underbrace{\sigma_m \cdot \eta_t}_{\text{stochastic innovation}}$$

where $a_{h,\ell}$ are state variables (lagged inflows) and $\eta_t$ is the sampled noise realization.

Whether the lagged inflows enter the Benders cuts is a separate, per-stage choice from whether they drive this inflow equation. A stage’s state_variables selection defaults to storage-only, so by default a PAR($p$) study’s cuts span storage alone — the inflow-lag dimensions are projected out, not zero-padded. Enabling inflow_lags on a stage restores the full storage-plus-lag cut; disabling it on every stage while fitting a PAR($p > 0$) model raises a model-quality warning, since the fitted inflow memory then never reaches any cut. The storage-only projection is the mechanism by which a PAR($p$) study can be coupled to a storage-only downstream boundary. See cut management.

4. Model Order Selection

The PAR order $p$ can vary by season. Cobre exposes one order-selection configuration key — estimation.order_selection — with exactly two accepted values: "pacf" (the classical periodic PACF method, the default) and "pacf_annual" (the same PACF rule augmented with the annual component, PAR(p)-A — see section 9). The classical method is described below. The fitting run is bounded by max_order (the largest lag considered, default 6), min_observations_per_season (groups with fewer observations are skipped), and an optional max_coefficient_magnitude safety net (any season whose fitted coefficient exceeds it is reduced to order 0 before the contribution analysis runs).

4.1 PACF (Periodic Partial Autocorrelation Function) — Default

The default method computes the periodic PACF via progressive periodic Yule-Walker matrix solves at orders $k = 1, 2, \ldots, p_{max}$, then selects the order using a significance threshold.

Algorithm:

1. For each order $k$ from 1 to $p_{max}$, build and solve the periodic Yule-Walker system (section 5.4) at order $k$. The last coefficient $\hat{\psi}^*_{m,k}$ from the order-$k$ solution is the periodic PACF value at lag $k$.

2. Select the order as the maximum lag with significant PACF:

   $$p_m = \max \left\{ k : |\text{PACF}_m(k)| > \frac{z_\alpha}{\sqrt{N_m}} \right\}$$

   where $z_\alpha = 1.96$ (95% confidence) and $N_m$ is the number of observations for season $m$. If no lag is significant, $p_m = 0$ (white noise).

3. Estimate AR coefficients at the selected order using the periodic Yule-Walker system (section 5.4).

Post-selection validation — iterative order reduction (Maceira & Damázio, 2006): After PACF selection, the recursively-composed contributions of each lag through the periodic monthly chain are computed. A negative composed contribution flags potential model instability — under SDDP the corresponding Benders cut can carry the negative composition into the future-cost recursion. When any season’s composed contribution is negative, the offending season’s AR ceiling is reduced and the PACF selection plus Yule-Walker fit are re-run at the new ceiling. The reduction iterates across all seasons until every season’s contribution recursion yields non-negative entries, or every offending season has been reduced to order 0.

For the PAR(p)-A path (section 9), two additional rules extend the PACF gate:

• Structural-zero short-circuit at lag 1. When the conditional FACP value at lag 1 is exactly zero — which happens when the standardised annual noise series collapses, typically because a degenerate HistoryClass::Constant or HistoryClass::Saturated bucket has zeroed the seasonal std (section 5.7) — the selected order is forced to 0 (white noise). This blocks degenerate buckets from injecting spurious AR structure.
• Minimum order 1 when lag 1 is well defined. When the lag-1 conditional FACP is non-zero but no lag exceeds the significance threshold, the model defaults to order 1 rather than order 0. Hydrological persistence makes a strict order-0 fit a poor default unless the lag-1 value is structurally absent.

5. Fitting Procedure

For multi-resolution studies (monthly→quarterly aggregation), the same fitting procedure applies after duration-weighted aggregation; see Multi-resolution studies.

This section documents the four-step procedure for fitting PAR(p) parameters from historical inflow data. The fitting is performed when the system derives parameters from inflow_history.parquet. When pre-computed parameters are provided directly in inflow_seasonal_stats.parquet and inflow_ar_coefficients.parquet, this procedure is not executed. The innovation scale is deliberately not part of the fitting: it is derived from the stored coefficients at load (section 3), identically for fitted and user-supplied models (section 5.6).

5.1 Notation

Let $Y_m = \{a_{h,t} : m(t) = m\}$ be the historical observations for season $m$. Define:

Symbol	Description
$N_m$	Number of observations for season $m$
$\bar{a}_m$	Sample mean for season $m$
$s_m$	Sample standard deviation for season $m$
$\gamma_m(\ell)$	Autocovariance at lag $\ell$ for season $m$
$\rho_m(\ell)$	Autocorrelation at lag $\ell$ for season $m$

5.2 Step 1 — Seasonal Means and Standard Deviations

Seasonal Mean:

$$\hat{\mu}_m = \bar{a}_m = \frac{1}{N_m} \sum_{t: m(t) = m} a_{h,t}$$

Seasonal Standard Deviation:

$$\hat{s}_m = \sqrt{\frac{1}{N_m} \sum_{t: m(t) = m} (a_{h,t} - \bar{a}_m)^2}$$

The estimator uses the population divisor $1/N_m$, not the Bessel-corrected $1/(N_m - 1)$. This matches the Maceira & Damázio (2006) convention and is shared by the classical PAR(p) and PAR(p)-A paths. The population divisor is required for self-consistent conditional FACP values and selected orders on the PAR(p)-A path — under a Bessel correction the sample-vs-population scale factor leaks through every Z⊗A cross-correlation. Using the same divisor for the classical path keeps the two paths’ seasonal-stats output reusable across configurations.

5.3 Step 2 — Seasonal Autocorrelations

The autocorrelation at lag $\ell$ for season $m$ is computed from standardized deviations.

Cross-seasonal autocovariance:

For observations at season $m$ with lag $\ell$ reaching back to season $m - \ell$ (mod $M$, where $M$ is the cycle length):

$$\hat{\gamma}_m(\ell) = \frac{1}{N_m^{(\ell)}} \sum_{t: m(t) = m} \left( a_{h,t} - \bar{a}_m \right) \left( a_{h,t-\ell} - \bar{a}_{m-\ell} \right)$$

where $N_m^{(\ell)}$ is the number of year-aligned valid pairs at lag $\ell$ for reference season $m$. The estimator uses the population divisor $1/N_m^{(\ell)}$, matching the convention adopted in section 5.2 and shared by the classical and PAR(p)-A paths.

Autocorrelation:

$$\hat{\rho}_m(\ell) = \frac{\hat{\gamma}_m(\ell)}{\hat{s}_m \cdot \hat{s}_{m-\ell}}$$

where $\hat{s}_{m-\ell}$ is the standard deviation of season $m - \ell$ (cyclically, so season 0 = season $M$).

5.4 Step 3 — Yule-Walker Equations

For each season $m$, the PAR(p) coefficients $\psi_{m,1}^*, \ldots, \psi_{m,p}^*$ in standardized form are found by solving the periodic Yule-Walker system. Unlike the classical (stationary) Yule-Walker equations where all rows use the same reference season, the periodic variant shifts the reference season per row. This correctly accounts for the non-Toeplitz covariance structure of periodic autoregressive processes.

Matrix construction: Entry $(i, j)$ (0-indexed, $0 \leq i,j < p$) is the correlation between the lagged observations $\tilde{a}_{t-(i+1)}$ and $\tilde{a}_{t-(j+1)}$; its reference season is that of the more recent of the two lags:

$$[\mathbf{R}_m]_{i,j} = \hat{\rho}_{\left(m - 1 - \min(i,j)\right) \bmod M}\bigl(|j - i|\bigr)$$

where $M$ is the number of seasons in the periodic cycle (e.g., 12 for monthly). The diagonal entries are always 1 (since $\hat{\rho}_{m'}(0) = 1$ for any season $m'$). The matrix is symmetric but not Toeplitz when $M > 1$, because the reference season varies across the entries.

RHS construction: The RHS holds the correlations of the target $\tilde{a}_t$ (season $m$) with each lag:

$$[\boldsymbol{r}_m]_i = \hat{\rho}_{m}(i + 1)$$

Equivalently: form the extended $(p{+}1) \times (p{+}1)$ correlation matrix of $(\tilde{a}_t, \tilde{a}_{t-1}, \ldots, \tilde{a}_{t-p})$; the system matrix is its lag block (rows/columns $1..p$) and the RHS is its first row (the target-lag correlations). Row $i$ of the system is the model’s second-moment recursion at lag $i{+}1$.

The full system is:

$$\begin{pmatrix} 1 & \hat{\rho}_{(m-1)}(1) & \hat{\rho}_{(m-1)}(2) & \cdots & \hat{\rho}_{(m-1)}(p{-}1) \\ \hat{\rho}_{(m-1)}(1) & 1 & \hat{\rho}_{(m-2)}(1) & \cdots & \hat{\rho}_{(m-2)}(p{-}2) \\ \hat{\rho}_{(m-1)}(2) & \hat{\rho}_{(m-2)}(1) & 1 & \cdots & \hat{\rho}_{(m-3)}(p{-}3) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \hat{\rho}_{(m-1)}(p{-}1) & \hat{\rho}_{(m-2)}(p{-}2) & \hat{\rho}_{(m-3)}(p{-}3) & \cdots & 1 \end{pmatrix} \begin{pmatrix} \psi_{m,1}^* \\ \psi_{m,2}^* \\ \psi_{m,3}^* \\ \vdots \\ \psi_{m,p}^* \end{pmatrix} = \begin{pmatrix} \hat{\rho}_{m}(1) \\ \hat{\rho}_{m}(2) \\ \hat{\rho}_{m}(3) \\ \vdots \\ \hat{\rho}_{m}(p) \end{pmatrix}$$

where all season indices are taken modulo $M$.

In matrix notation: $\mathbf{R}_m \boldsymbol{\psi}_m^* = \boldsymbol{r}_m$

where:

• $\mathbf{R}_m$ is the $p \times p$ periodic correlation matrix (symmetric but not Toeplitz for $M > 1$)
• $\boldsymbol{r}_m$ is the vector of target autocorrelations with per-row reference season shifting

Note

For a single-season model ($M = 1$), all rows use the same reference season and the matrix reduces to the classical Toeplitz Yule-Walker matrix. The periodic formulation is the general case that correctly handles multi-season (e.g., monthly) data.

Solution:

$$\hat{\boldsymbol{\psi}}_m^* = \mathbf{R}_m^{-1} \boldsymbol{r}_m$$

The system is solved via Gaussian elimination with partial pivoting (for small systems with $p \leq 10$, this is numerically adequate).

Reference formulation

This is the classical periodic Yule-Walker system of the PAR literature (Hipel & McLeod, 1994, ch. 14), written in correlation form. The canonical presentation states it in autocovariances $\gamma$ with original-unit coefficients $\psi_{m,j}$; dividing each row $k$ of that system by $s_m\,s_{m-k}$ and substituting $\psi^*_{m,j} = \psi_{m,j}\,s_{m-j}/s_m$ yields exactly the matrix above. The lag-0 row of the same derivation — the variance identity $\gamma_0^{(m)} = \sum_j \psi_{m,j}\,\gamma_j^{(m)} + \sigma_m^2$ — divided by $s_m^2$ becomes the unit-variance decomposition $1 = \sum_j \psi^*_{m,j}\,\rho_m(j) + r_m^2$ that section 3’s closure evaluates on the model’s implied ACF.

5.5 Step 4 — Store Standardized Coefficients

The Yule-Walker solution $\psi_{m,\ell}^*$ is in standardized form — the direct output of step 3. It is stored as-is in inflow_ar_coefficients.parquet. No conversion to original units is performed, and no other quantity is written alongside it: the innovation scale is not a fitting output. It is derived from $\psi^*$ afterward, uniformly for every model regardless of its origin, by the periodic-ACF closure of section 3 — not computed here.

When every season of the cycle shares the same AR order, that closure value coincides exactly with this step’s own Yule-Walker system:

$$r_m = \sqrt{1 - \boldsymbol{\psi}_m^{*\top} \boldsymbol{r}_m}$$

(using this section’s $\boldsymbol{\psi}_m^*$ and the system’s RHS vector $\boldsymbol{r}_m$ of section 5.4). When per-season orders differ across the cycle, the two need not agree: the closure solves for the model’s implied ACF jointly across every season (section 3), not from this season’s sample autocorrelations alone.

Only $\psi^*_{m,\ell}$ (one row per lag) is written to inflow_ar_coefficients.parquet.

5.6 Innovation Scale — Derived at Load, Not Fitted

The innovation standard deviation for season $m$ is recovered at load from the closure-derived innovation scale $r_m$ (see section 3):

$$\hat{\sigma}_m = \hat{s}_m \cdot r_m$$

When every season of the cycle shares the same AR order, section 5.5 showed this reduces to the fitting system’s own closed-form expression:

$$\hat{\sigma}_m = \hat{s}_m \sqrt{1 - \boldsymbol{r}_m^\top \mathbf{R}_m^{-1} \boldsymbol{r}_m}$$

5.7 Historical Bucket Classification

Before the seasonal stats and AR coefficients are used by the order-selection rules, each per-(hydro, season) historical bucket is classified by the shape of its observations. The classification can override the empirical $(\hat{\mu}_m, \hat{s}_m)$ for fitting purposes, and the override propagates to both the classical PAR(p) and the PAR(p)-A paths because both paths share the seasonal-stats producer.

Four classes are defined:

Class	Detection rule	Override applied
Default	None of the conditions below	None — use empirical $(\hat{\mu}_m, \hat{s}_m)$
Constant	Every observation equals the same value within float tolerance	$(\hat{\mu}_m, \hat{s}_m) \leftarrow (\text{value}, 0)$
ManyNegative	Strictly negative observations exceed 10% of the bucket	None — diagnostic only, fit proceeds on the empirical stats
Saturated	The modal value (rounded to m³/s) occupies more than 50% of observations	$(\hat{\mu}_m, \hat{s}_m) \leftarrow (\text{cap}, 0)$

The classifier runs in the priority order Constant → ManyNegative → Saturated → Default. Constancy takes precedence over negative-pathology detection, which in turn takes precedence over saturation.

Why a zero seasonal std short-circuits the fit

When the override sets $\hat{s}_m = 0$ for a season, every downstream fitter degenerates predictably:

• On the classical PAR(p) path, the periodic autocorrelation $\hat{\rho}_m(\ell)$ becomes zero by the zero-std guard in section 5.3, so the PACF selection (section 4.1) reports no significant lag and returns order 0 implicitly.
• On the PAR(p)-A path, the standardised noise series collapses, the conditional FACP at lag 1 evaluates to exactly zero, and the structural-zero short-circuit (section 4.1) returns order 0 explicitly.

Either way, the bucket cannot inject spurious autoregressive structure into adjacent months’ PACFs, and no spatial-correlation contribution flows from it during scenario generation.

Interpretation of each class

• Constant captures plants whose incremental inflow is structurally constant for a given month — typically regulated or transposed flows where the upstream subtraction yields the same value every year. Forcing $(\text{value}, 0)$ records the deterministic level without inventing autoregressive dynamics.
• Saturated captures flow caps (turbine or reservoir capacity) and low-flow constants (transposed ecological flows). The modal value is treated as the cap. There is no magnitude threshold — a cap of 0 m³/s qualifies just as readily as a cap at installed capacity.
• ManyNegative flags buckets that the upstream incremental-inflow construction has driven below zero for more than 10% of observations. The condition is recorded for operator diagnostics but does not override the fit — the cause is upstream-data quality, not a methodological signal.
• Default is the standard path; the empirical stats and the chosen order-selection rule decide the order.

5.8 Partial-year studies and the pre-study lag window

A study horizon may be narrower than the seasonal cycle — e.g. a monthly model ($M = 12$) running only September–December. The per-season fitting described above must then handle seasons that have few or no in-window observations. Two rules keep it well-defined.

Lag-reachability. A season is lag-reachable only if some stage of the (extended) horizon carries it. Each historical observation is resolved to a season from the stage date ranges, falling back to the season-map calendar for dates predating the horizon; an observation whose resolved season has no stage at all is skipped — its statistics would never be consumed. Full-cycle history therefore does not perturb a partial-year fit.

Pre-study lag synthesis (for $p > 0$). The first study stage’s autoregressive lags reach back to seasons before the study start. For each lag $k = 1, \ldots, \min(p,\, M - 1)$, the season $k$ calendar positions before the first study season is introduced as a pre-study season (modular on the true cycle length $M$) — unless that season is already covered by a study stage (an in-window wrap lag, handled by the cycle-correct lag lookup). The seasonal statistics $(\hat{\mu}_m, \hat{s}_m)$ of those out-of-window seasons are estimated from history exactly as for in-window seasons, then feed the lag terms of the opening study stages — both the coefficient conversion $\psi_{m,\ell} = \psi^*_{m,\ell}\, \hat{s}_m / \hat{s}_{m-\ell}$ and the deterministic base.

The wrap uses the true cycle length $M$ (the number of seasons in the season map), not the number of seasons in the study window, together with a season offset equal to the season of the first study stage — so, e.g., a March-start study maps the lag-1 season to February, not December.

Full-cycle invariance. When the study spans the full cycle (every season already has a study stage) or carries no out-of-window history, nothing is synthesized and the fit is bit-identical to before.

6. Validation Invariants

Cobre enforces a small, well-defined set of invariants; they fall into three groups by where the check runs. Stability is a fitting-time guard rather than a post-hoc root test.

Enforced when loading inflow_ar_coefficients.parquet:

1. AR order derivation: the number of coefficient rows per (hydro_id, stage_id) determines the AR order $p_m$, and the lags must be present and contiguous $\{1, 2, \ldots, p_m\}$.
2. Stationarity of directly-supplied coefficients: when the AR coefficients for a (hydro, season) group are supplied directly rather than produced by the internal fitting procedure (section 5), the periodic-ACF closure (section 3) must accept them — the closure system must be solvable (non-singular), every implied innovation variance $r_m^2 = 1 - \sum_\ell \psi^*_{m,\ell}\,\rho_m(\ell)$ must be strictly positive, every implied autocorrelation must satisfy $|\rho_m(k)| \le 1$, and the periodic monodromy — the product, taken once around the full seasonal cycle, of each season’s AR companion matrix — must have spectral radius strictly below 1 (assessed via a conservative upper-bound estimate, so borderline sets are rejected rather than accepted). A set failing any condition is rejected outright. A model whose coefficients come from the internal fitting procedure is exempt from this check: its own order-reduction procedure (section 4.1) already guards stability during the fit.

Enforced when validating the assembled model:

3. Positive sample std (error): a season with AR order $> 0$ must have $s_m > 0$ — a zero seasonal std cannot normalise the AR coefficients.
4. High explained variance (warning): $R^2 = 1 - r_m^2 > 0.99$ (the AR part explains more than 99% of the season’s variance — a likely over-fit) raises an advisory model-quality warning; it is not fatal.

Enforced during fitting, not as a post-hoc test:

5. Stability: rather than testing the roots of the AR polynomial, the fitter guards stability through the Maceira & Damázio contribution-based order reduction (section 4.1) — a negative lag-1 coefficient, a negative composed contribution, or a coefficient exceeding max_coefficient_magnitude reduces the season’s order and refits. The periodic Yule-Walker solve additionally truncates the order when its correlation matrix is singular (section 5.4), so a rank-deficient system never yields coefficients.

7. PAR-to-LP Transformation

This section derives the explicit algebraic transformation from the canonical PAR(p) model (section 1) into the form consumed by the LP subproblem. The derivation identifies three precomputable components that are cached once at initialization and reused at every forward-pass stage transition.

7.1 Canonical Standardized Form

The PAR(p) model (section 1) operates on deviations from the seasonal mean, normalised by the seasonal (marginal) standard deviation $s_{m(t)}$ — the same basis the stored coefficients $\psi^*$ use, and the same normalisation as the reference formulation (standardise by the marginal std, not the innovation std). In this standardized form:

$$\frac{a_{h,t} - \mu_{m(t)}}{s_{m(t)}} = \sum_{\ell=1}^{p} \psi^*_{m(t),\ell}\, \frac{a_{h,t-\ell} - \mu_{m(t-\ell)}}{s_{m(t-\ell)}} + r_{m(t)}\,\varepsilon_t$$

where:

• $\psi^*_{m(t),\ell}$: the stored AR coefficients, standardized by the seasonal std $s_m$ (the direct Yule-Walker output of section 5.4)
• $s_{m(t)}$: the seasonal (marginal) standard deviation for season $m(t)$
• $r_{m(t)} = \sigma_{m(t)}/s_{m(t)} \in (0, 1]$: the innovation standard deviation of the standardized process — the closure-derived innovation scale of section 3, $\sqrt{1 - \sum_\ell \psi^*_{m,\ell}\,\rho_{m(t)}(\ell)}$, with $\rho_{m(t)}$ the model’s implied periodic ACF
• $\varepsilon_t \sim \mathcal{N}(0, 1)$: unit-variance innovation noise

The innovation of the $s_m$-standardized process is not unit-variance: because that process has unit marginal variance, its one-step innovation has standard deviation $r_{m(t)} \le 1$ — which is exactly why $r_{m(t)}$ appears here explicitly. The next step converts $\psi^*_{m,\ell}$ to original-unit $\psi_{m,\ell}$ for use in the LP.

7.2 Coefficient Conversion

The stored standardized coefficients $\psi^*_{m,\ell}$ are converted to original-unit coefficients $\psi_{m,\ell}$ at runtime using the seasonal standard deviations from inflow_seasonal_stats.parquet:

$$\psi_{m,\ell} = \psi^*_{m,\ell} \cdot \frac{s_m}{s_{m-\ell}}$$

The innovation standard deviation is also derived at this preprocessing step, from the closure-derived innovation scale $r_m$ (section 3):

$$\sigma_m = s_m \cdot r_m$$

These conversions are performed once at LP construction time. They require only the seasonal stats ($s_m$), the stored $\psi^*_{m,\ell}$, and the closure-derived $r_m$ — no historical data.

7.3 LP-Ready Form

Multiplying both sides of the canonical form (section 7.1) by $s_{m(t)}$ and rearranging yields the LP-ready equation (the noise term becomes $s_{m(t)}\,r_{m(t)}\,\varepsilon_t = \sigma_{m(t)}\,\varepsilon_t$):

$$a_{h,t} = \sum_{\ell=1}^{p} \psi_{m(t),\ell} \cdot a_{h,t-\ell} + \left[ \mu_{m(t)} - \sum_{\ell=1}^{p} \psi_{m(t),\ell} \cdot \mu_{m(t-\ell)} \right] + \sigma_{m(t)} \cdot \varepsilon_t$$

where $\psi_{m(t),\ell}$ and $\sigma_{m(t)}$ are derived from stored quantities as described in section 7.2.

This decomposes the inflow into three additive components:

1. Lag contribution: $\displaystyle\sum_{\ell=1}^{p} \psi_{m(t),\ell} \cdot a_{h,t-\ell}$ — linear function of past inflows (state variables or known values)
2. Deterministic base: $\displaystyle\mu_{m(t)} - \sum_{\ell=1}^{p} \psi_{m(t),\ell} \cdot \mu_{m(t-\ell)}$ — constant offset per (stage, hydro), precomputed once
3. Stochastic innovation: $\sigma_{m(t)} \cdot \varepsilon_t$ — noise draw scaled by the seasonal innovation standard deviation

7.4 Deterministic Base

The deterministic base is defined as:

$$b_{h,m(t)} = \mu_{m(t)} - \sum_{\ell=1}^{p} \psi_{m(t),\ell} \cdot \mu_{m(t-\ell)}$$

This is a precomputed constant per (stage, hydro) pair. It absorbs the mean-adjustment arithmetic that would otherwise be repeated at every forward-pass stage transition. With this definition, the LP-ready form (section 7.3) simplifies to:

$$a_{h,t} = \sum_{\ell=1}^{p} \psi_{m(t),\ell} \cdot a_{h,t-\ell} + b_{h,m(t)} + \sigma_{m(t)} \cdot \varepsilon_t$$

For partial-year studies, the lag-season means $\mu_{m(t-\ell)}$ for seasons preceding the study start are sourced from the pre-study lag window (section 5.8); when no such statistic exists for a given lag season, that lag’s mean contribution is treated as zero.

7.5 LP RHS Patching Operation

The lagged inflows $a_{h,t-\ell}$ are LP variables, not substituted values. In the LP (see LP Formulation), they appear with coefficients $-\psi_{m(t),\ell}$ in the AR dynamics constraint row, and separate equality constraints fix each lag variable to its incoming state value:

$$a_{h,t-\ell} = \hat{a}_{h,t-\ell}$$

where $\hat{a}_{h,t-\ell}$ is patched per scenario to the actual lagged inflow from the trajectory record.

Because the lag contribution $\sum_\ell \psi \cdot a_{h,t-\ell}$ is carried by the constraint matrix (not the RHS), the AR dynamics constraint RHS reduces to:

$$\text{RHS}_{h,t} = b_{h,m(t)} + \sigma_{m(t)} \cdot \varepsilon_t$$

where:

• $b_{h,m(t)}$ is the deterministic base for (stage, hydro), precomputed once at LP construction (section 7.4)
• $\sigma_{m(t)}$ is the noise scale for (stage, hydro), derived from the closure-derived $r_m$ (section 3) at initialization (section 7.2)
• $\varepsilon_t$ is the scenario noise draw for this (stage, hydro)

The $\psi_{m(t),\ell}$ coefficients are written into the constraint matrix once at LP construction time as the coefficients on the lagged inflow variables; they are not recomputed per scenario.

No division, no mean subtraction, no repeated coefficient transformation — the three precomputed LP components eliminate all redundant arithmetic from the hot path.

7.6 Summary of LP Components

Component	Symbol	Shape per stage	LP Role	Source
Lag coefficients	$\psi_{m(t),\ell}$	One per (hydro, lag)	Constraint matrix (AR dynamics row)	Derived from stored $\psi^*$ and $s_m$ at initialization (section 7.2)
Deterministic base	$b_{h,m(t)}$	One per hydro	AR dynamics constraint RHS (fixed term)	Precomputed from $\mu$ and $\psi$
Noise scale	$\sigma_{m(t)}$	One per hydro	AR dynamics constraint RHS (noise factor)	Derived from the closure-derived $r_m$ (section 3) and $s_m$ at initialization (section 7.2)

8. Spatial Correlation Factorisation

The PAR(p) fitting procedure (section 5) produces per-hydro noise terms $\varepsilon_t$ that are treated as independent across hydro plants. Generating spatially correlated scenarios requires factorising the cross-hydro correlation matrix $C$ so that a vector of independent standard normal draws can be mapped to correlated noise. This section documents the choice of factorisation method and the rationale.

8.1 The Problem with Cholesky

The classical approach applies Cholesky factorisation: given $C = L L^\top$ with $L$ lower-triangular, correlated noise is obtained as $L z$ where $z \sim \mathcal{N}(0, I)$. Cholesky requires $C$ to be strictly positive-definite. In practice, estimated correlation matrices from hydro inflow series are frequently near-singular or rank-deficient for two reasons:

• Short sample records: Brazilian hydro inflow series span only a few decades of record, yielding a historical record length $N$ that is comparable to the number of hydro plants in some subsystems. When $N$ is close to the matrix dimension, the sample eigenvalues of $C$ cluster near zero.
• Heterogeneous series: Plants with near-identical hydrological regimes (upstream–downstream pairs, same river basin) produce columns that are nearly linearly dependent, reducing the effective rank of $C$ below its nominal dimension.

A near-singular $C$ causes Cholesky to fail or to produce numerically degenerate lower triangular factors. A separate filtering pass to remove “degenerate” hydros would be required before the factorisation, discarding information and introducing a non-transparent pre-processing decision.

8.2 Eigendecomposition with Clipped Square Root

Cobre uses the symmetric matrix square root via eigendecomposition. The correlation matrix is decomposed as:

$$C = U \Lambda U^\top$$

where $U$ is the orthogonal matrix of eigenvectors and $\Lambda = \mathrm{diag}(\lambda_1, \ldots, \lambda_n)$ is the diagonal matrix of eigenvalues. The symmetric square root is then:

$$C^{1/2} = U \Lambda^{1/2} U^\top$$

To handle near-singular matrices, any eigenvalue $\lambda_i < 0$ (arising from floating-point rounding in the sample estimate) is clipped to zero before taking the square root:

$$\tilde{\Lambda}^{1/2} = \mathrm{diag}\!\left(\sqrt{\max(\lambda_1, 0)},\, \ldots,\, \sqrt{\max(\lambda_n, 0)}\right)$$

Clipping negative eigenvalues to zero is the spectral projection of the sample matrix onto the positive-semidefinite cone — the nearest positive-semidefinite matrix in Frobenius norm; see Higham (2002) for the related nearest-correlation-matrix problem.

Correlated noise is then generated as $C^{1/2} z$ where $z \sim \mathcal{N}(0, I)$.

8.3 Why Eigendecomposition

The spectral form handles rank-deficient correlation matrices natively: eigenvectors corresponding to clipped (zero) eigenvalues contribute nothing to the factorisation, which is the correct behaviour for directions of zero variance. No prior filtering of degenerate hydro plants is needed.

The clipping threshold acts as a single, transparent parameter controlling which near-zero eigenvalues are treated as structural zeros. The cross-entity correlation structure is preserved for all eigenvalues above the threshold.

8.4 Trade-offs

Property	Eigendecomposition (Cobre)	Cholesky
Handles rank-deficient $C$	Yes — clipping makes it robust	No — requires positive-definiteness
Computational cost	Higher (full eigendecomposition)	Lower on well-conditioned matrices
Degenerate-hydro filtering pass	Not required	Required for near-singular $C$
Transparency of approximation	Single clipping threshold	Opaque numerical failure or pivot

The higher computational cost is acceptable because the factorisation is performed once per study configuration and not on the hot path of the forward pass.

9. Annual Component Extension (PAR(p)-A)

The classical PAR(p) of section 1 captures temporal dependence at lags up to a small order $p$ (kept low for monthly cycles, since the periodic Yule-Walker system becomes ill-conditioned at higher orders). On long Brazilian hydro series this is enough to reproduce the within-year persistence but not the multi-year persistence visible in dry/wet super-periods of the historical record. The PAR(p)-A extension adds a single annual coefficient on top of the periodic AR structure to capture that longer-range persistence without inflating the AR order.

The extension is selected by the order-selection method pacf_annual. When active, the model carries one additional triple per (hydro, season) on top of the classical parameter set.

9.1 Augmented Model

Let $A_{h,t-1}$ denote the rolling 12-month average of incremental inflows ending one stage before $t$:

$$A_{h,t-1} = \frac{1}{12} \sum_{j=1}^{12} a_{h,\, t-j}$$

The PAR(p)-A model augments section 1 with the standardised deviation of $A_{h,t-1}$ from its own seasonal mean:

$$a_{h,t} \;=\; \mu_{m(t)} \;+\; \sum_{\ell=1}^{p} \psi_{m(t),\ell}\,(a_{h,t-\ell} - \mu_{m(t-\ell)}) \;+\; \hat{\psi}_{m(t)}\,(A_{h,t-1} - \mu^A_{m(t)}) \;+\; \sigma_{m(t)} \cdot \varepsilon_t$$

where:

• $\mu^A_{m}$, $\sigma^A_{m}$: sample mean and population-divisor standard deviation of the season’s own annual regressor — the values of $A_{h,t-1}$ over stages $t$ in season $m$, i.e. the rolling windows whose most recent observation falls in the preceding season (section 9.3)
• $\hat{\psi}_{m(t)}$: original-unit annual coefficient at season $m(t)$ — derived at runtime from the standardised stored coefficient (section 9.4)
• All other symbols carry their classical meaning from section 1

When the PAR(p)-A extension is inactive, the annual term is absent and the model reduces exactly to section 1.

9.2 Annual Component Parameters

For each (hydro, season) the PAR(p)-A path stores three additional quantities:

Quantity	Symbol	Description
Standardised annual coefficient	$\psi$	Yule-Walker output for the annual term — dimensionless
Annual seasonal mean	$\mu^A_m$	Sample mean of season $m$‘s annual regressor $A_{h,t-1}$ (m³/s)
Annual seasonal std	$\sigma^A_m$	Population-divisor std of season $m$‘s annual regressor $A_{h,t-1}$ (m³/s, $> 0$)

The standardised coefficient $\psi$ is the direct output of the extended periodic Yule-Walker system below (section 9.5). Storage of $\mu^A_m$ and $\sigma^A_m$ alongside the seasonal statistics of $a_{h, \cdot}$ enables the runtime unit conversion of section 9.4 without re-reading the historical record.

9.3 Estimating the Annual Seasonal Statistics

Form every rolling 12-month average $A_t = \frac{1}{12} \sum_{j=0}^{11} a_{h,\, t - 11 + j}$ the chronological history admits, and assign each window to the season following its most recent observation — the season whose stages use that window as their regressor $A_{h,t-1}$. The statistics stored for season $m$ therefore describe exactly the regressor that season’s stages see. For each (hydro, season $m$) bucket of values $\{A^{(i)}\}$:

$$\hat{\mu}^A_m \;=\; \frac{1}{N^A_m} \sum_{i} A^{(i)} \qquad \hat{\sigma}^A_m \;=\; \sqrt{\frac{1}{N^A_m} \sum_{i} \bigl(A^{(i)} - \hat{\mu}^A_m\bigr)^2}$$

Both estimators use the population divisor $1/N^A_m$, matching the convention of section 5.2 and ensuring no sample-vs-population scale factor leaks into the conditional FACP of section 9.5. At least 13 chronological observations are required for a hydro to participate in PAR(p)-A — that is the minimum needed to form one rolling 12-month average.

9.4 Runtime Unit Conversion

The stored standardised coefficient $\psi$ is converted to the original-unit coefficient $\hat{\psi}$ at LP construction time using the seasonal stats and annual stats:

$$\hat{\psi}_{m} \;=\; \psi \cdot \frac{s_m}{\sigma^A_m}$$

The conversion mirrors section 7.2 for the classical AR coefficients. The annual term $\hat{\psi}_{m(t)} \cdot \bigl(A_{h,t-1} - \mu^A_{m(t)}\bigr)$ is then expanded through the window definition at LP construction: each of the 12 lag coefficients gains $\hat{\psi}_{m(t)}/12$ on top of its classical value, and the deterministic base absorbs the corresponding mean contribution through the lag-season means — the LP carries no separate annual variable, and the rolling-window value is realized entirely by the lag state variables it already carries.

9.5 Order Selection and Coefficient Estimation

PAR(p)-A order selection conditions on the annual noise series. The order-selection input is the conditional FACP at lag $k$, defined as the partial autocorrelation between the standardised current-season residual and the standardised residual at lag $k$, conditioned on the intermediate standardised annual noise series $Z$ and the previous annual innovation $A_{t-1}$. Computing the conditional FACP requires a partitioned covariance decomposition that distinguishes $Z \otimes Z$, $Z \otimes A$, and $A \otimes Z_{-1}$ blocks.

The conditional FACP feeds the PACF order-selection rule of section 4.1, with the two PAR(p)-A-specific extensions (structural-zero short-circuit and minimum-order-1) already described there.

Cross-covariance divisor

The $Z \otimes Z$ block uses the same year-aligned population divisor as the classical autocovariance (section 5.3). The $Z \otimes A$ and $A \otimes Z_{-1}$ blocks use a max-bucket-size divisor:

$$\hat{\gamma}_{Z \otimes A}(\ell) \;=\; \frac{1}{\max(|A|,\, |Z|)} \sum_i \bigl(Z^{(i)} - \bar{Z}\bigr)\bigl(A^{(i)} - \bar{A}\bigr)$$

The max-bucket convention is required because $A$ excludes the first year of $Z$ by construction (a rolling 12-month window cannot anchor in the first 12 observations). The strict-pair count would distort the scale of the cross-correlations and bias the conditional FACP. The PAR(p) path never uses Z⊗A cross-correlations, so the divisor question is PAR(p)-A-specific.

Extended periodic Yule-Walker

Once the order $p$ is selected, the coefficients $(\psi^*_{m,1}, \ldots, \psi^*_{m,p}, \psi)$ are recovered by solving the extended periodic Yule-Walker system:

$$\mathbf{R}^{\,\text{ext}}_m \begin{pmatrix} \psi^*_{m,1} \\ \vdots \\ \psi^*_{m,p} \\ \psi \end{pmatrix} = \boldsymbol{r}^{\,\text{ext}}_m$$

where $\mathbf{R}^{\,\text{ext}}_m$ is the $(p+1) \times (p+1)$ partitioned covariance whose first $p$ rows replicate the classical periodic Yule-Walker rows (section 5.4) and whose last row adds the $Z \otimes A$ and $A \otimes A$ entries. The RHS $\boldsymbol{r}^{\,\text{ext}}_m$ appends the $A \otimes Z_{-1}$ target.

As in the classical case, the innovation scale is not a fitting output. For a PAR(p)-A model the closure of section 3 runs on the effective 12-lag system: the annual regressor is a linear functional of the last 12 inflows, so expanding it per lag yields an effective periodic AR whose standardized lag-$\tau$ coefficient is the classical $\psi^*_{m,\tau}$ (zero beyond $p_m$) plus the annual contribution $\psi \cdot s_{m-\tau} / (12\,\sigma^A_m)$. The innovation scale $r_m$ derives from that effective system’s implied ACF, and the runtime noise scale is $\sigma_m = s_m \cdot r_m$ exactly as in section 3.

9.6 Iterative Order Reduction

The Maceira & Damázio iterative reduction of section 4.1 is applied across the full periodic cycle on the PAR(p)-A path as well: after the initial fit, the recursively-composed contributions of each AR lag through the periodic monthly chain are evaluated, and any negative-contribution season has its AR ceiling reduced before refit. The annual coefficient $\psi$ does not enter the contribution-chain check — it is anchored to the rolling annual mean and so does not propagate through the lag chain.

9.7 Activation and Fallback

Configuration	Path
order_selection: "pacf"	Classical PAR(p) — annual triple absent (section 5)
order_selection: "pacf_annual"	PAR(p)-A — annual triple required for every (hydro, season)
Bucket flagged HistoryClass::Constant or Saturated (section 5.7)	Effective order 0 on either path; annual term suppressed when the seasonal std collapses
Hydro with fewer than 13 observations on PAR(p)-A path	Hard failure during fitting (no silent fallback to classical)

The two paths share the seasonal-stats producer of section 5.2; switching between them does not silently change $\hat{\mu}_m$ or $\hat{s}_m$. The PAR(p)-A path uses the same spatial-correlation factorisation as the classical path (section 8); the extension affects only the temporal model.

9.8 Conditioning-Swap Exactness

Section 3 shows that re-conditioning the seasonal stats ($\mu_m$, $s_m$) alone — leaving $\psi^*$ and the closure-derived $r_m$ untouched — rescales the classical dynamics while preserving the correlation structure exactly, for any per-season rescaling. The PAR(p)-A annual term does not share this property in general.

Like $r_m$, the annual seasonal std $\sigma^A_m$ is a functional of the process’s own second-moment structure — it is not an independent parameter of the annual term. Unlike $r_m$, it is carried as stored conditioning data (section 9.2) rather than re-derived from $\psi^*$ at load. The runtime annual coefficient $\hat\psi_m$ (section 9.4) is $\hat\psi_m = \psi \cdot s_m / \sigma^A_m$, so a conditioning swap preserves the annual term’s contribution to the correlation structure only when $\sigma^A_m$ rescales by the same factor as $s_m$ in every season — a uniform rescaling $\tilde s_m = \lambda\, s_m$ for all $m$, giving $\tilde\sigma^A_m = \lambda\,\sigma^A_m$ and leaving $\hat\psi_m$ unchanged. A season-varying conditioning swap — rescaling $s_m$ by a different factor in different seasons — leaves $\sigma^A_m$ at its pre-swap value while $s_m$ moves, shifting $\hat\psi_m$ and perturbing the preserved correlation structure.

Implementation in Cobre

The methodology above defines the PAR(p) inflow model; the tabs below cover how Cobre’s software surface configures, feeds, and reports on it.

Configure

Cobre’s config.json top-level estimation block, the scenarios/ PAR(p) input files, and scenarios/correlation.json configure the estimation and correlation pipeline the methodology above describes. This tab is the field-level configuration reference; the equations these fields feed are in the sections above.

config.json — estimation Block

The estimation block is a top-level key of config.json (a sibling of training), and controls the PAR(p) fitting pipeline whenever inflow_history.parquet is present:

{
  "estimation": {
    "max_order": 6,
    "order_selection": "pacf",
    "min_observations_per_season": 30,
    "max_coefficient_magnitude": null
  }
}

Field	Type	Default	Description
max_order	integer	6	Maximum lag order p_max considered by PACF order selection (methodology §4.1).
order_selection	string	"pacf"	Order-selection criterion. Exactly two values are accepted: "pacf" (classical PAR(p), methodology §1–§8) and "pacf_annual" (PAR(p)-A, methodology §9).
min_observations_per_season	integer	30	Recommended minimum observation count per (hydro, season) group. Below this count, estimation still proceeds — see Implementation notes.
max_coefficient_magnitude	float or null	null	Optional safety net on fitted coefficient magnitude — see Implementation notes for the exact effect.

Setting "order_selection": "pacf_annual" activates the PAR(p)-A annual component (methodology §9): the periodic Yule-Walker system is extended with the annual cross-correlation term, per-season sample statistics of the rolling 12-month average are computed, and the fitted triple (ψ, annual mean, annual std) is written to inflow_annual_component.parquet — see the Inputs & Outputs tab. When absent (or set to "pacf"), only the classical PAR(p) path runs and no annual triple is produced.

Inflow Source Resolution — Choosing Which Files to Provide

The PAR(p) inflow model is built from up to five files in scenarios/. Three of them drive path resolution — their presence or absence selects which of seven estimation paths Cobre executes:

Symbol	File	Role
H	scenarios/inflow_history.parquet	Raw historical observations for fitting
S	scenarios/inflow_seasonal_stats.parquet	User-supplied seasonal mean/std
R	scenarios/inflow_ar_coefficients.parquet	User-supplied AR coefficients

The other two files layer orthogonally on top of the resolved path: scenarios/correlation.json wins on every path when present (identity correlation otherwise, unless a path estimates it from residuals); and scenarios/inflow_annual_component.parquet is honored only on the pass-through paths (rows 2, 3, 7 below) — the estimation paths (4, 5, 6) always overwrite it with fitted values.

#	H	S	R	Path	Seasonal stats	AR coefficients	Correlation
1	0	0	0	Deterministic	no PAR model	none	identity, unless correlation.json provided
2	0	1	0	UserStatsWhiteNoise	user file	order-0 (white noise)	identity, unless correlation.json provided
3	0	1	1	UserProvidedNoHistory	user file	user file	identity, unless correlation.json provided
4	1	0	0	FullEstimation	fitted from H	fitted from H (PACF + Yule-Walker + Maceira & Damázio)	estimated from H residuals, unless correlation.json provided
5	1	0	1	UserArHistoryStats	fitted from H	user file	estimated from H residuals using user coefficients, unless correlation.json provided
6	1	1	0	PartialEstimation	user file (used only for the Yule-Walker solve)	fitted from H	estimated from H residuals using the user stats, unless correlation.json provided
7	1	1	1	UserProvidedAll	user file	user file	identity, unless correlation.json provided (history is not re-consumed)

Cases with R = 1 but H = 0 and S = 0 collapse to row 1 — AR coefficients alone cannot drive estimation.

Practical recipes:

Goal	Files to provide	Path
Smoke-test the LP without stochasticity	(no scenarios files)	1
Deterministic seasonal levels, no autoregression	inflow_seasonal_stats.parquet	2
Fully user-specified PAR(p) without raw observations	inflow_seasonal_stats.parquet, inflow_ar_coefficients.parquet	3
Hands-off: fit everything from raw observations	inflow_history.parquet	4
Fit stats from history, override the AR structure	inflow_history.parquet, inflow_ar_coefficients.parquet	5
Override the levels (mean/std) but let Cobre fit the AR	inflow_history.parquet, inflow_seasonal_stats.parquet	6
Provide every parameter, including the PAR(p)-A annual term	All three of H, S, R (and optionally the annual file)	7
Pin a custom spatial correlation on any path	Add correlation.json	any

scenarios/inflow_history.parquet — Raw Historical Observations

Required to drive any of the estimation paths (H = 1 above). One row per historical observation:

Column	Type	Required	Description
hydro_id	INT32	Yes	Hydro plant identifier (matches id in hydros.json).
date	DATE	Yes	Observation date.
value_m3s	DOUBLE	Yes	Historical inflow observation (m³/s); must be finite.

Every hydro plant in hydros.json must have at least one observation in this file when estimation is active. When a season has fewer observations than min_observations_per_season, estimation still proceeds — see Implementation notes for the exact behavior.

scenarios/correlation.json — Spatial Correlation Configuration

Named correlation profiles, each grouping entities with a symmetric correlation matrix, plus an optional stage-to-profile schedule:

{
  "method": "spectral",
  "profiles": {
    "default": {
      "correlation_groups": [
        {
          "name": "basin_south",
          "entities": [
            { "type": "inflow", "id": 0 },
            { "type": "inflow", "id": 1 }
          ],
          "matrix": [
            [1.0, 0.7],
            [0.7, 1.0]
          ]
        }
      ]
    },
    "wet_season": {
      "correlation_groups": [
        {
          "name": "basin_south",
          "entities": [
            { "type": "inflow", "id": 0 },
            { "type": "inflow", "id": 1 }
          ],
          "matrix": [
            [1.0, 0.85],
            [0.85, 1.0]
          ]
        }
      ]
    }
  },
  "schedule": [{ "stage_id": 0, "profile_name": "wet_season" }]
}

Field	Required	Description
method	Yes	Correlation decomposition method. "spectral" (eigendecomposition + clipped square root, methodology §8.2) is preferred; "cholesky" is also accepted for existing case files and is treated identically.
profiles	Yes	Named profiles, each holding one or more correlation_groups. Must not be empty.
schedule	No	Optional array of { stage_id, profile_name } entries activating a named profile for a specific stage. Absent means the "default" profile (or the sole profile) applies to every stage. Every profile_name referenced must exist in profiles.

Each correlation group’s matrix must be square, match the length of entities, have diagonal entries exactly 1.0, off-diagonal entries in [-1.0, 1.0], and be symmetric.

Valid entity types in entities[].type:

• "inflow" — hydro inflow series (id matches id in hydros.json)
• "load" — stochastic load demand (id matches id in buses.json)
• "ncs" — non-controllable source availability (id matches id in non_controllable_sources.json)

Same-type enforcement: every entity within a single correlation group must share the same type. Mixing types in one group is rejected at case-load time; to correlate inflow with load, define separate groups (one per type) with matching correlation structure.

Inputs & Outputs

This is a topic-scoped index of the files PAR(p) estimation touches — it names each file and its role, it does not repeat their field-by-field schemas beyond what the Configure tab already shows for the config-authoring surface. The exhaustive, field-by-field case-directory and output reference is owned by the Reference corpus (Case Directory Format and Output Format pages).

Inputs

File	Role
scenarios/inflow_history.parquet	Raw historical inflow observations per (hydro, date) — drives the estimation path (H) when present (see the Configure tab).
scenarios/inflow_seasonal_stats.parquet	User-supplied or fitted seasonal statistics per (hydro, stage). The std_m3s column is the stored sample standard deviation of the historical series — not the innovation standard deviation used for noise scaling (see Implementation notes).
scenarios/inflow_ar_coefficients.parquet	User-supplied or fitted standardized AR coefficients, one row per (hydro, stage, lag) — schema [hydro_id, stage_id, lag, coefficient]. The AR order is derived from the row count, not stored explicitly.
scenarios/inflow_annual_component.parquet	Optional user-supplied PAR(p)-A annual triple (coefficient, mean, std of the rolling 12-month average). Honored only on the pass-through paths — the estimation paths overwrite it with fitted values (see the Configure tab).
scenarios/correlation.json	Named spatial-correlation profiles and optional stage schedule. Wins over estimated or identity correlation on every inflow-source path when present.
initial_conditions.json (hydro past_inflows rows)	The AR-lag seed: past inflow values ordered most-recent-first, consumed by the opening study stages’ lag terms (methodology §7.5).

For the complete field-by-field schema of each file above, see the Case Directory Format reference page in the Reference corpus.

Externally-fitted PAR models

Cobre accepts fully user-supplied parameters — you may fit a PAR($p$) model with any method outside Cobre and hand it the two parquet files (the UserProvidedNoHistory and UserProvidedAll paths on the Configure tab). Cobre then skips fitting, validates the inputs, and derives the innovation scale $r_m$ from your supplied $\psi^*$ via the periodic-ACF closure (see Deriving the innovation scale). The fields do not encode Cobre’s fitting method — they encode the two-plane representation — so any fit maps onto them by the table below. Let $s_m$ be your series’ seasonal sample std and $\psi_{m,\ell}$ your original-unit AR coefficients.

Column (file)	Fill with	Meaning
mean_m3s (inflow_seasonal_stats)	$\mu_m$	Seasonal mean of the series (m³/s)
std_m3s (inflow_seasonal_stats)	$s_m$	Seasonal sample std of the series (m³/s)
coefficient (inflow_ar_coefficients)	$\psi^*_{m,\ell} = \psi_{m,\ell}\,s_{m-\ell}/s_m$	AR coefficient standardized by the seasonal std

There is no column for your model’s own innovation std $\sigma_m$: Cobre always derives $r_m = \sigma_m/s_m$ from the $\psi^*_{m,\ell}$ you supply, the same way it derives $r_m$ for an internally-fitted model (see Deriving the innovation scale). If your own fit’s $\hat\sigma_m$ is not exactly consistent with your supplied $\psi^*$ under the periodic-ACF closure, Cobre’s derived $\sigma_m$ differs from your original fit’s — the closure derives the noise scale so the seasonal std $s_m$ you supply is honored exactly.

Rules Cobre enforces (or that silently mis-scale the model if broken):

• Standardize by the seasonal std, consistently. The coefficients are those of the series standardized by its seasonal sample std — the classical standardization — and the $s_m$ you divide by must be exactly the std_m3s you store: Cobre reconstructs $\psi$ and $\sigma$ from the stored std_m3s, so a mismatch is silent.
• Cobre hard-rejects a coefficient set that fails its stationarity gate. The gate runs the periodic-ACF closure on your supplied $\psi^*$ and rejects the whole load when the closure system is singular, the implied innovation variance is not strictly positive for any season, any implied autocorrelation magnitude exceeds 1, or a conservative spectral-radius estimate of the periodic monodromy reaches 1. The gate runs only on directly-supplied coefficients — coefficients Cobre fits internally keep the fitter’s own fallbacks.
• Lags must be contiguous $1, 2, \ldots, p$; the AR order $p$ is implicit in the row count (there is no order column). White noise = no coefficient rows for that (hydro, stage), i.e. order 0, for which the closure gives $r_m = 1$.
• Resolvable season context is required whenever any AR order is $>0$. Every (hydro, stage) group with at least one coefficient row must resolve to a season via season_definitions or a per-stage season_id (see Case Directory Format and Running Cobre — configuration) — this applies on every load path, including a plain file-load with no fitting.

Outputs

File	Role
output/stochastic/inflow_seasonal_stats.parquet	Fitted seasonal statistics, same schema as the input file — written whenever the estimation path ran.
output/stochastic/inflow_ar_coefficients.parquet	Fitted AR coefficients, same schema as the input file — written whenever the estimation path ran.
output/stochastic/inflow_annual_component.parquet	Fitted PAR(p)-A annual triple — written only when order_selection: "pacf_annual" and the estimation path ran.
output/stochastic/fitting_report.json	Per-hydro summary of the selected AR order and coefficients — a human-readable complement to the parquet outputs above.
output/stochastic/correlation.json	Round-trip export of the resolved correlation model, in the same schema as the scenarios/correlation.json input.
output/stochastic/noise_openings.parquet	The backward-pass opening tree actually used during training, exported when exports.stochastic: true. Copying it back into scenarios/noise_openings.parquet replays the same opening tree on a subsequent run — see Scenario Generation for the opening-tree mechanics themselves (out of scope here).

For the complete output schema (columns, types, file layout), see the Output Format reference page in the Reference corpus.

Implementation notes

Non-normative software behavior for PAR(p) estimation — what Cobre does at runtime, beyond the equations above. This tab references the methodology body for the derivations rather than restating them.

std_m3s vs sigma_m

The two files in the Inputs & Outputs tab store two different standard deviations, and the software surface never conflates them: std_m3s (column in inflow_seasonal_stats.parquet) is the stored sample standard deviation of the historical series for a (hydro, season) — a fixed, on-disk number. sigma_m, the innovation standard deviation that scales the stochastic innovation term, is never stored; it is derived at load from std_m3s and the innovation scale r_m, which Cobre always derives from the AR coefficients via the periodic-ACF closure (methodology §3). Reading std_m3s as if it were sigma_m — or vice versa — silently rescales every noise draw.

Where the closure runs relative to the fit

The closure derivation of r_m is a load step, not a fitting step: it runs after model assembly on every input path — pure file-load and all estimation paths alike — so a fitted model and a user-supplied model with the same coefficients get the same innovation scale. When per-season AR orders differ across the periodic cycle, the closure-derived value can differ slightly from the single-season Yule-Walker identity evaluated on the sample autocorrelations (methodology §5.5) — the closure solves for the model’s own implied ACF jointly across all seasons, which is what makes the generated seasonal marginal std match std_m3s exactly. When every season shares the same order, the two coincide to double precision.

Model-quality diagnostic: explained variance

The model-quality warning for a near-deterministic season reports an explained_variance field ($R^2 = 1 - r_m^2$), firing when $R^2 > 0.99$ (methodology §6). It is a library-surface diagnostic type only: Cobre computes it internally and does not render it on the CLI or expose it through the python bindings. Keep this $R^2$ distinct from r_squared, the unexplained variance reported in the stationarity gate’s rejection message — the two are complementary quantities describing different situations: explained_variance is a quality signal on an accepted model, while r_squared in a rejection message is the reason a model was rejected outright.

Order-selection behavior

The classical PACF order-selection rule (methodology §4.1) and the Maceira & Damázio iterative reduction it triggers are applied automatically whenever Cobre estimates from inflow_history.parquet. Two runtime details sit alongside that rule but are not part of the derivation:

• min_observations_per_season is a warning, not a load-time error. When a (hydro, season) group has fewer historical observations than the configured value (default 30), Cobre logs a model-quality warning naming the plant and season and proceeds with the fit anyway. The threshold is a quality signal, not a hard gate.
• max_coefficient_magnitude is an independent safety net. When set and any fitted coefficient for a (hydro, season) exceeds the threshold in magnitude, Cobre clears that group’s coefficients entirely (order reset to 0, for which the closure then derives r_m = 1 — pure white noise). This runs as a separate pass from the Maceira & Damázio contribution-chain reduction (methodology §4.1) — a coefficient can trigger this magnitude cutoff even when its composed contribution is non-negative.

History classification runs before any order selection

Before PACF selection sees a (hydro, season) bucket, the history classification pipeline (methodology §5.7) runs automatically: constant, saturated, and negative-dominated series are detected and, for the first two, routed to a degenerate order-0 fit rather than left to over-fit a structurally uninformative bucket. This classification is unconditional — it always runs when estimating from history, on both the classical PAR(p) and PAR(p)-A paths.

PAR(p)-A annual-component runtime behavior

Setting order_selection: "pacf_annual" (Configure tab) extends the estimation pipeline with the steps of methodology §9. At runtime this adds one requirement not visible in the equations: a hydro needs at least 13 chronological observations to participate in PAR(p)-A — the minimum to form one rolling 12-month average (§9.3). A hydro with fewer observations on the PAR(p)-A path is a hard estimation failure; there is no silent fallback to the classical path (methodology §9.7).

Determinism guarantee

PAR(p) fitting is one of the model families covered by Cobre’s order-stable parallel fitting guarantee: coefficients are fit independently per hydro and reassembled into a canonical per-hydro slot, so the fitted model is a function of the inputs alone — independent of how many threads run or how hydros are distributed across them. See Determinism Guarantees §3 for the full mechanism, which this chapter’s estimation pipeline shares with the computed-FPHA fitting pipeline described in Hydro Production Models.

Cross-References

• LP Formulation — AR inflow dynamics in the LP: state expansion, lag column pinning, reduced-cost extraction
• Inflow Non-Negativity — Methods for handling negative realizations produced by the PAR(p) model
• Scenario Generation — When external scenarios are used in training, a PAR model is fitted to the external data for backward pass opening tree generation. The fitting procedure (section 5 above) applies equally to this derived model.
• Notation Conventions — Defines inflow symbols ($a_{h,t}$, $\mu_m$, $\psi_{m,\ell}$, $\sigma_m$) and unit conventions