PAR(p) Inflow Model

Purpose

This chapter defines the Periodic Autoregressive model of order $p$ (PAR(p)) used to capture temporal correlation in inflow time series. It states the model, its parameters and the quantities stored for it (§1); the form the LP stage subproblem consumes (§2); the estimation of the parameters from historical inflow data, including order selection (§3); the closure that derives the innovation scale from the AR coefficients (§4); the validation invariants (§5); the estimation and factorisation of the spatial correlation across hydros (§6); and the optional PAR(p)-A extension, which adds a single annual coefficient on top of the periodic AR structure to capture multi-year hydrological persistence (§7).

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 §4)
• $\varepsilon_t$ $\sim \mathcal{N}(0, 1)$: Innovation (standardized noise)
• $m(t)$: The resolved season of stage $t$; lag seasons are calendar predecessors (see Notation Conventions)

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

1.1 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 (§4)

1.2 Stored and Derived 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 (§4.1), 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 through the coefficient conversion of §2.2 and the innovation standard deviation of §4.2, 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 quantities. For each hydro and season the data model stores:

Stored quantity	Symbol	Description
Seasonal sample mean	$\mu_m = \bar{a}_m$	Mean of historical observations for season $m$
Seasonal sample std	$s_m$	Standard deviation of historical observations for season $m$
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 number of stored coefficients of each (hydro, stage).

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

Relation to other PAR(p) implementations

Other PAR(p) implementations may follow different conventions:

• they often write the autoregressive coefficients as $\phi$, which this chapter writes $\psi$ (original units) and $\psi^*$ (standardized);
• some store the innovation standard deviation, while this model stores the seasonal sample standard deviation $s_m$ and derives the innovation standard deviation $\sigma_m = s_m\, r_m$ at load (§4);
• some draw residuals from a lognormal distribution, while this model’s innovations are Gaussian and negative realizations are treated as in Inflow Non-Negativity Solution Methods.

For equivalent terms in other planning tools, see the Glossary.

2. LP Form

This section derives the explicit algebraic transformation from the canonical PAR(p) model (§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.

2.1 Canonical Standardized Form

The PAR(p) model (§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 §3.5)
• $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 §4.1, $\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.

2.2 Coefficient Conversion

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

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

The innovation standard deviation $\sigma_m$ of §4.2 is also derived at this preprocessing step, from the closure-derived innovation scale $r_m$ (§4.1).

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.

2.3 LP-Ready Form

Multiplying both sides of the canonical form (§2.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 §2.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

2.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 (§2.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 (§3.8); when no such statistic exists for a given lag season, that lag’s mean contribution is treated as zero.

2.5 Lagged Inflows Pinned by Column Bounds

The lagged inflows $a_{h,\ell}$ are LP variables, not substituted values. In the LP, they appear with coefficients $-\psi_{m(t),\ell}$ in the AR dynamics constraint row, and each incoming lag column is pinned to the trajectory’s lag value $\hat{a}_{h,\ell}$ by setting both of its column bounds to that value: there are no fixing rows, as for every state variable in State Augmentation §2. The value $\hat{a}_{h,\ell}$ is patched per scenario from the trajectory record.

Because the lag contribution $\sum_\ell \psi_{m(t),\ell} \cdot a_{h,\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 (§2.4)
• $\sigma_{m(t)}$ is the noise scale for (stage, hydro), derived from the closure-derived $r_m$ (§4.1) at initialization (§2.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.

Whether the lagged inflows enter a stage’s cuts is a per-stage choice; a cut that projects them out carries their contribution in its intercept at the trial lag values (see State Augmentation §7).

2.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 (§2.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$ (§4.1) and $s_m$ at initialization (§2.2)

3. Estimation

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

This section documents the procedure for fitting PAR(p) parameters from historical inflow data. The fitting is performed when the system derives parameters from the inflow history. When pre-computed seasonal statistics and AR coefficients are provided directly, this procedure is not executed. The innovation scale is deliberately not part of the fitting: it is derived from the stored coefficients at load (§4.1), identically for fitted and user-supplied models (§4.2).

The figure below traces both origins of the coefficients to the runtime model. When the coefficients are fitted from the inflow history, the seasonal statistics of each (hydro, season) are estimated and the bucket classification applied (§3.2, §3.3); a season whose seasonal std is zero takes order 0 directly. Every other season selects its order from the periodic PACF (§3.4, §3.6) and solves the periodic Yule-Walker system at that order for its standardized coefficients $\psi^*$ (§3.5). A season is reset to order 0 when its first coefficient is negative, after the initial fit or any re-fit, or when its initial fit has a coefficient of magnitude above the optional magnitude bound. A season with a negative composed influence loses one lag of its order ceiling and is re-selected and re-fitted, until no season fails; the coefficients are then stored (§3.7), and the residual correlation is estimated when the case supplies none (§6.1). Supplied coefficients pass the stationarity gate of invariant 2 (§5) instead, and a set that fails it is rejected. When the case supplies only the coefficients, or only the seasonal statistics, alongside the history, the history provides the other part and the residual correlation is estimated as above. Every model, fitted or supplied, then passes through the closure that derives its innovation scale $r_m$ (§4.1), which rejects a singular closure system or a negative implied residual variance; the runtime coefficients $\psi$ and noise scale $\sigma$ follow (§4.2).

3.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$

3.2 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.

3.3 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)$
Many negative	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 → many negative → 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, the fit short-circuits that season on both paths:

• On the classical PAR(p) path, the season takes order 0 directly; its PACF (§3.6) is not evaluated.
• On the PAR(p)-A path, the season likewise takes order 0 directly, with annual coefficient $\psi^{A*}_m = 0$; its conditional FACP (§7.5) is not evaluated.

The seasonal statistics that feed the fit are estimated from the history, even when the case supplies its own, and estimating them needs at least two observations per season, so no season reaches the fit with a nonzero std and $N_m < 2$. A season with a single observation ($N_m = 1$) rejects the case. A season with none ($N_m = 0$) has no statistics; the fit reads its std as zero and short-circuits it in the same way on both paths. On the PAR(p)-A path it also does so when the season’s annual-regressor bucket is empty ($N^A_m = 0$) or has zero std ($\hat{\sigma}^A_m = 0$, §7.3), and the season then carries no annual contribution.

The zero-std guard of §3.4 also zeroes every sample autocorrelation that involves the season, so 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.
• Many negative 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.

3.4 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 §3.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$).

Zero-std guard. $\hat{\rho}_m(\ell) = 0$ when $\hat{s}_m$ or $\hat{s}_{m-\ell}$ is zero, or when season $m$ has no year-aligned pair at lag $\ell$ ($N_m^{(\ell)} = 0$). Every estimate is clamped to $[-1, 1]$.

3.5 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)$, $1 \leq i, j \leq p$, is the correlation between the lagged observations $\tilde{a}_{t-i}$ and $\tilde{a}_{t-j}$; its reference season is that of the more recent of the two lags:

$$[\mathbf{R}_m]_{i,j} = \hat{\rho}_{\left(m - \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:

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

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$.

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 season indices are taken modulo $M$ on $\{1, \ldots, M\}$, so season $0$ is season $M$ (see Notation Conventions).

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

where:

• $\mathbf{R}_m$ is the $p \times p$ periodic correlation matrix (symmetric but not Toeplitz for $M > 1$)
• $\hat{\boldsymbol{\rho}}_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} \hat{\boldsymbol{\rho}}_m$$

The system is solved via Gaussian elimination with partial pivoting, which is numerically adequate for the small systems that periodic Yule-Walker estimation solves.

Singular system. When the periodic Yule-Walker system is singular at order $k$, the PACF of §3.6 stops at lag $k-1$, so only lags below $k$ can be selected and the fit at the selected order meets no singular system.

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,\ell}$; dividing each row $\ell'$ of that system by $s_m\,s_{m-\ell'}$ and substituting $\psi^*_{m,\ell} = \psi_{m,\ell}\,s_{m-\ell}/s_m$ yields exactly the matrix above. The lag-0 row of the same derivation — the variance identity $\gamma_m(0) = \sum_\ell \psi_{m,\ell}\,\gamma_m(\ell) + \sigma_m^2$ — divided by $s_m^2$ becomes the unit-variance decomposition $1 = \sum_\ell \psi^*_{m,\ell}\,\rho_m(\ell) + r_m^2$ that the closure of §4.1 evaluates on the model’s implied ACF.

3.6 Order Selection

The PAR order $p$ can vary by season. Two order-selection rules are available: the classical periodic PACF rule, described below, and the order-selection rule with the annual component, which augments the same PACF rule for PAR(p)-A (§7). The fit considers lags up to a maximum lag $p_{max}$ and applies an optional magnitude bound to the coefficients of the initial fit (below).

The classical rule computes the periodic PACF (periodic partial autocorrelation function) 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 (§3.5) 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_{0.975}}{\sqrt{N_m}} \right\}$$

   where $z_{0.975} = 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 (§3.5).

The figure below applies the rule to a synthetic series: a twelve-season PAR(1) process with unit marginal variance in every season, simulated with a fixed random seed for $N_m = 60$ years and shown for season 6. The upper panel compares the sample periodic autocorrelation $\hat{\rho}_6(\ell)$ of §3.4 with the model’s implied autocorrelation of §4.1, which for a PAR(1) is $\rho_m(\ell) = \prod_{j=0}^{\ell-1} \psi^*_{m-j,1}$ with season indices taken modulo $M$. The lower panel shows the periodic PACF of step 1, for the sample and for the model, against the band $\pm z_{0.975}/\sqrt{N_m}$, which is the significance threshold of step 2: a lag is significant when its sample PACF lies outside the band. The model’s PACF is $\psi^*_{m,1}$ at lag 1 and zero beyond, so the autocorrelation decays over all six lags while the PACF cuts off after lag 1. Here the sample PACF at lag 1 lies outside the band and the sample PACF at lags 2 to 6 lies inside it, so, with $p_{max} = 6$, the rule selects $p_6 = 1$.

Post-selection validation — iterative order reduction (Maceira & Damázio, 2006): After PACF selection, the composed influence of each lag on the current inflow through the periodic chain is computed for every season, holding the inflows before $t-\ell$ and the innovations fixed:

$$\frac{\partial a_{h,t}}{\partial a_{h,t-\ell}} = \psi_{m(t),\ell} + \sum_{\ell'=1}^{\ell-1} \psi_{m(t),\ell'}\, \frac{\partial a_{h,t-\ell'}}{\partial a_{h,t-\ell}}, \qquad \ell = 1, \ldots, p_{m(t)}$$

where $\psi_{m,\ell}$ is the original-unit coefficient of §2.2, taken as zero when $s_{m-\ell} = 0$ or $\ell > p_m$, and each derivative on the right follows from the same recursion at stage $t-\ell'$. The check is a sign condition on the composed lag influence, not a stationarity test: a season fails when any of its composed influences is negative. A failing season loses one lag of its order ceiling (initially $p_{max}$) and is re-selected and re-fitted at the new ceiling; a season whose ceiling reaches 0 takes order 0, and the loop repeats until no season fails. A season whose first coefficient satisfies $\psi^*_{m,1} < 0$, after the initial fit or after any re-fit, is reset to order 0, and so is a season whose initial fit has a coefficient of magnitude above the optional magnitude bound.

For the PAR(p)-A path (§7), 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 inflow series $Z$ of §7.5 collapses in the preceding season, typically because a degenerate constant or saturated bucket has zeroed that season’s std (§3.3) — 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.

3.7 Stored Standardized Coefficients

The Yule-Walker solution $\psi_{m,\ell}^*$ is in standardized form — the direct output of §3.5. It is stored as-is, one coefficient per lag. 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 §4.1 — not computed here.

3.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 $\ell = 1, \ldots, \min(p,\, M - 1)$, the season $\ell$ calendar positions before the first study season is introduced as a pre-study season — 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 (§2.2) and the deterministic base (§2.4).

The walk steps back through the calendar occurrences of the first study season’s own resolution, whatever the seasons’ numbering — 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 equals the full-cycle fit.

4. Innovation-Scale Closure

4.1 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}$ (§3.5) to every lag $\ell \ge 1$, not just $\ell \le p_m$:

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

Solving this system jointly across every season $m = 1, \ldots, M$ and every lag $\ell = 1, \ldots, P_h$ (where $P_h = \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 (§3.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.

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

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

(using the Yule-Walker solution $\boldsymbol{\psi}_m^*$ and the RHS vector $\hat{\boldsymbol{\rho}}_m$ of §3.5). 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, not from this season’s sample autocorrelations alone.

4.2 Runtime Quantities

From the stored quantities and the derived $r_m$, the LP requires two additional quantities computed once at initialization: the original-unit AR coefficients $\psi_{m,\ell}$ of §2.2 (for LP constraint matrix entries) and the innovation standard deviation (the noise scale), recovered at load from the closure-derived innovation scale $r_m$ of §4.1:

$$\sigma_m = s_m \cdot r_m$$

No further autocorrelation values are needed beyond the closure derivation of §4.1. 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 coefficients, with no historical data and no separate noise-scale input.

For a model estimated from history this is $\hat{\sigma}_m = \hat{s}_m\, r_m$, and when every season of the cycle shares the same AR order, the equal-order identity of §4.1 reduces it to the fitting system’s own closed-form expression:

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

5. Validation Invariants

Novomodelo enforces a small, well-defined set of invariants; they fall into three groups by where the check runs. Fitted coefficients are not tested for stationarity; supplied coefficients pass the stationarity gate of invariant 2.

Enforced when loading the AR coefficients:

1. AR order derivation: the number of stored coefficients per (hydro, stage) 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 (§3), the periodic-ACF closure (§4.1) 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 lie above a numerical floor, every implied autocorrelation must satisfy $|\rho_m(\ell)| \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. Independent of origin, however, the closure derivation itself (§4.1) is a final load-time guard: a singular closure system is a hard load error naming the hydro, and a non-finite derived $r_m$ — a negative implied residual variance $1 - \sum_\ell \psi^*_{m,\ell}\,\rho_m(\ell)$ from a degenerate near-unit-root fit — is a hard load error naming the hydro and season.

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.

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

4. Non-negative composed influence: after order selection, every season’s composed lag influences (§3.6) must be non-negative; a failing season loses one lag of its order ceiling and is re-selected and re-fitted until no season fails. A season whose first coefficient satisfies $\psi^*_{m,1} < 0$, after the initial fit or after any re-fit, is reset to order 0, and so is a season whose initial fit has a coefficient of magnitude above the optional magnitude bound. On the classical path, a singular periodic Yule-Walker system bounds the selectable order below its own order (§3.5), so the classical fit never solves a singular system.

6. Spatial Correlation Factorisation

The residuals of the fitted autoregression are correlated across hydro plants: §6.1 estimates that correlation, and the rest of the section factorises it. 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.

6.1 Correlation Estimation

When the inflow model is fitted from the inflow history, in whole or in part, and the case supplies no correlation, the correlation matrix is estimated from the historical residuals of the autoregression. The standardized residual of hydro $h$ at historical period $t$ is

$$\tilde{\varepsilon}_{h,t} = \tilde{a}_{h,t} - \sum_{\ell=1}^{p_{m(t)}} \psi^*_{m(t),\ell}\, \tilde{a}_{h,t-\ell}$$

where $\tilde{a}_{h,t} = (a_{h,t} - \hat{\mu}_{m(t)}) / \hat{s}_{m(t)}$ is the inflow standardized by the seasonal statistics of §3.2, set to zero when $\hat{s}_{m(t)} = 0$, and $\psi^*_{m,\ell}$ are the model’s standardized coefficients (§1.2). A period contributes a residual only when all $p_{m(t)}$ of its lagged inflows are observed. On the PAR(p)-A path the residual uses the lag coefficients only: the annual term of §7.1 does not enter it.

Two estimates are formed from these residuals: the pooled estimate $\hat{C}$ from the residuals of every season, and the per-season estimate $\hat{C}_m$ from the residuals of the periods in season $m$. Each off-diagonal entry is the Pearson correlation of the two hydros’ residuals over the periods at which both hydros have a residual (a pairwise-complete estimate), with the covariance and the standard deviations taken with the Bessel divisor, one less than the number of such periods, unlike the population divisor of §3.2. The entry is zero when fewer than two such periods exist or when either residual standard deviation is below machine epsilon, and it is clamped to $[-1, 1]$; every diagonal entry is one.

The stages of season $m$ use $\hat{C}_m$ when every pair of hydros has at least 30 such periods in season $m$; all other stages use $\hat{C}$. Because each entry is computed over its own set of periods, neither estimate need be positive semidefinite; §6.3 handles an indefinite matrix.

6.2 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: hydro inflow records often span only a few decades, yielding a historical record length $N^{\text{hist}}$ that is comparable to the number of hydro plants in some subsystems. When $N^{\text{hist}}$ 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.

6.3 Eigendecomposition with Clipped Square Root

Novomodelo 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$ (which the pairwise-complete estimate of §6.1 can produce) 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)$.

With the clipped root in $C^{1/2}$, the generated noise has covariance $\tilde{C} = U \max(\Lambda, 0)\, U^\top$, where $\max(\Lambda, 0) = \mathrm{diag}(\max(\lambda_1, 0), \ldots, \max(\lambda_n, 0))$; this is the projection above, and it equals $C$ when no eigenvalue is negative. Because $C$ has a unit diagonal, $\sum_i U_{hi}^2 \lambda_i = 1$ for every hydro $h$, with $U_{hi}$ the component of eigenvector $i$ on hydro $h$, so

$$\tilde{C}_{hh} = \sum_i U_{hi}^2 \max(\lambda_i, 0) = 1 + \sum_{i:\, \lambda_i < 0} U_{hi}^2 \lvert \lambda_i \rvert \ge 1$$

No renormalisation follows, so clipping raises the marginal variance of each affected hydro’s noise above one.

6.4 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 rule is fixed and transparent — exactly the strictly negative eigenvalues are set to zero, nothing else is altered. The eigen-directions with non-negative eigenvalues keep their eigenvalues unchanged. When any eigenvalue is clipped, the generated noise has covariance $\tilde{C}$ (§6.3) instead of $C$: its diagonal exceeds one for every hydro with a component along a clipped eigenvector, and its off-diagonal entries in general differ from those of $C$.

6.5 Trade-offs

Property	Eigendecomposition (Novomodelo)	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	Clips exactly the negative eigenvalues; the generated variance can exceed one	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.

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

The classical PAR(p) of §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 hydro inflow 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 (Treistman et al., 2020) 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 rule with the annual component (§3.6). When active, the model carries one additional triple per (hydro, season) on top of the classical parameter set.

7.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_{\ell=1}^{12} a_{h,\, t-\ell}$$

The window always averages the last twelve observations, so $A_{h,t-1}$ is a one-year average only on a monthly cycle ($M = 12$); with the extension active on any hydro, the AR dynamics row (§2.5) spans twelve lags, or the largest classical order if that is larger.

The PAR(p)-A model augments §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)}) \;+\; \psi^A_{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 (§7.3)
• $\psi^A_{m(t)}$: original-unit annual coefficient at season $m(t)$ — derived at runtime from the standardised stored coefficient (§7.4)
• All other symbols carry their classical meaning from §1

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

7.2 Annual Component Parameters

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

Quantity	Symbol	Description
Standardised annual coefficient	$\psi^{A*}_m$	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^{A*}_m$ is the direct output of the extended periodic Yule-Walker system below (§7.5). Storage of $\mu^A_m$ and $\sigma^A_m$ alongside the seasonal statistics of $a_{h, \cdot}$ enables the runtime unit conversion of §7.4 without re-reading the historical record.

7.3 Estimating the Annual Seasonal Statistics

Form every rolling 12-month average $A_t = \frac{1}{12} \sum_{\ell=1}^{12} a_{h,\, t + 1 - \ell}$ 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^{(w)}\}$:

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

Both estimators use the population divisor $1/N^A_m$, matching the convention of §3.2 and ensuring no sample-vs-population scale factor leaks into the conditional FACP of §7.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.

7.4 Runtime Unit Conversion

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

$$\psi^A_{m} \;=\; \psi^{A*}_m \cdot \frac{s_m}{\sigma^A_m}$$

The conversion mirrors §2.2 for the classical AR coefficients. The annual term $\psi^A_{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 $\psi^A_{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.

7.5 Order Selection and Coefficient Estimation

PAR(p)-A order selection conditions on the annual regressor. Let $Z$ be the standardised inflow series, the inflow standardised by its seasonal statistics (§3.2), and let the annual regressor $A_{t-1}$ of §7.1 be standardised by its own seasonal statistics (§7.3). The order-selection input is the conditional FACP at lag $k$, defined as the partial autocorrelation between the current-stage value of $Z$ and its value at lag $k$, conditioned on the intermediate lags of $Z$ and on $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 §3.6, 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 (§3.4). 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_w \bigl(Z^{(w)} - \bar{Z}\bigr)\bigl(A^{(w)} - \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^{A*}_m)$ 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^{A*}_m \end{pmatrix} = \hat{\boldsymbol{\rho}}^{\,\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 (§3.5) and whose last row adds the $Z \otimes A$ and $A \otimes A$ entries. The RHS $\hat{\boldsymbol{\rho}}^{\,\text{ext}}_m$ appends the $A \otimes Z_{-1}$ target.

Singular extended system. The extended matrix can be singular even when its leading $p \times p$ block, the classical periodic Yule-Walker matrix of §3.5, is not, so the guarantee of invariant 4 (§5) covers the classical path only. A season whose extended system is singular at the selected order takes order 0 with annual coefficient $\psi^{A*}_m = 0$.

As in the classical case, the innovation scale is not a fitting output. For a PAR(p)-A model the closure of §4.1 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-$\ell$ coefficient is the classical $\psi^*_{m,\ell}$ (zero beyond $p_m$) plus the annual contribution $\psi^{A*}_m \cdot s_{m-\ell} / (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 §4.2.

7.6 Iterative Order Reduction

The Maceira & Damázio iterative reduction of §3.6 is applied across the full periodic cycle on the PAR(p)-A path as well: after the initial fit, the composed influence of each AR lag through the periodic chain is evaluated, and any season with a negative composed influence has its AR ceiling reduced before refit. The annual coefficient $\psi^{A*}_m$ does not enter the composed-influence check — it is anchored to the rolling annual mean and so does not propagate through the lag chain.

7.7 Activation and Fallback

Condition	Path
Classical order-selection rule	Classical PAR(p) — annual triple absent (§3)
Order-selection rule with the annual component	PAR(p)-A — annual triple required for every (hydro, season)
Bucket classified constant or saturated (§3.3)	Effective order 0 on either path; annual term suppressed when the seasonal std collapses
Extended periodic Yule-Walker system singular at the selected order (§7.5)	PAR(p)-A season at order 0 with annual coefficient $\psi^{A*}_m = 0$
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 §3.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 (§6).

7.8 Conditioning-Swap Exactness

The two-plane split of §1.2 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 (§7.2) rather than re-derived from $\psi^*$ at load. The runtime annual coefficient $\psi^A_m$ (§7.4) is $\psi^A_m = \psi^{A*}_m \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 = c\, s_m$ for all $m$, giving $\tilde\sigma^A_m = c\,\sigma^A_m$ and leaving $\psi^A_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 $\psi^A_m$ and perturbing the preserved correlation structure.

Implementation in Novomodelo

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

Configure

Novomodelo’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 shows which PAR inputs to supply and how the estimation and correlation are configured; 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 on the estimation paths (rows 4–6 of the table in the next section):

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

The estimation section of the configuration reference lists every field and its default, and states the min_observations_per_season behaviour.

Setting "order_selection": "pacf_annual" activates the PAR(p)-A annual component (methodology §7): 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 ($\psi^{A*}_m$, 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 Novomodelo executes:

Symbol	File	Role
H	scenarios/inflow_history.parquet	Windowed 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 (the fit uses statistics estimated from H)	fitted from H	estimated from H residuals, 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.

An estimated correlation (rows 4–6, when no correlation.json is provided) holds a pooled "default" profile. In a cycle of more than one season it also holds one season_<id> profile for each season whose every hydro pair has at least 30 residual pairs — the season id zero-padded to the width of the largest id, as in season_03 for a twelve-season cycle — and a schedule that maps each stage of such a season to its profile. Every other stage uses the pooled profile (methodology §6.1).

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 a history file	inflow_seasonal_stats.parquet, inflow_ar_coefficients.parquet, and, for an AR order above 0, recent_observations covering the lag slots	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 Novomodelo 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 — Windowed Historical Observations

The windowed record is required to drive any of the estimation paths (H = 1 above) and, whenever present, is the base layer the derived inflow-lag seed casts from (see “Lag and accumulator seeding” below).

Every hydro plant in hydros.json must have at least one observation window in this file when estimation is active. The columns, window rules and roles of the file are in scenarios/inflow_history.parquet.

Lag and accumulator seeding — recent_observations

Independently of which estimation path (if any) inflow_history.parquet drives, its windowed record is the seed for the opening study stage’s PAR lag slots and mid-period inflow accumulator. initial_conditions.json’s optional recent_observations array conditions that seed: entries carry the same hydro_id / start_date / end_date / value_m3s shape as inflow_history.parquet and shadow it day-wise wherever the two overlap (the two series are never blended). value_m3s accepts negative values in both files: the quantity is incremental inflow, so a negative window is real hydrology (see Inflow Non-Negativity Solution Methods).

When scenarios/inflow_ar_coefficients.parquet is supplied, the record must fully cover each lag slot up to the maximum AR order and, with an annual component, each slot of the twelve-period annual window that the study’s own completed periods do not replace; otherwise the case is refused at load (BusinessRuleViolation). A gap in a deeper slot, which the study’s own completed periods push out of the lag window before the study ends, only draws a warning.

When the coefficients are estimated, this coverage check does not run: a lag period the record covers only in part seeds from the mean over its covered days, and a period with no covered day seeds to zero.

On either path, a first-stage period already in progress at the study start that the record covers only in part is accepted, with a warning; with no covered day, its accumulator seeds to zero.

scenarios/correlation.json — Spatial Correlation Configuration

Named profiles of correlation groups, with 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" }]
}

Correlation is modelled independently per entity class — inflow, load, and NCS each carry their own spatial correlation structure, and there is no cross-class correlation term (see Scenario Generation §2.1).

Every field and load rule of the file is in scenarios/correlation.json.

Out-of-sample fallback

Under the out_of_sample forward scheme, a stage whose sampling_method is historical_residuals or selective draws plain Monte Carlo noise, and the run logs a warning that names the entity class and the affected stage ids. For the message, see Scenario Generation — Implementation notes tab, Warnings.

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 Format and Output Format pages).

Inputs

File	Role
scenarios/inflow_history.parquet	Windowed historical inflow observations per (hydro, [start_date, end_date)). It has four roles, from estimation to the historical-residual openings (listed below the table).
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 (recent_observations rows)	Conditioning layer over the windowed inflow_history record: observed inflow for partial periods before the study start, shadowing inflow_history day-wise wherever the two overlap. Feeds the derived inflow-lag seed consumed by the opening study stages’ lag terms and mid-period accumulator (methodology §2.5, and the Configure tab).

scenarios/inflow_history.parquet has four roles:

1. It drives the estimation path (H) when present.
2. It is the base layer for the derived inflow-lag seed on every path (see the Configure tab).
3. It supplies the historical replay windows of the historical forward scheme.
4. It supplies the historical-residual openings of the stages whose sampling_method is historical_residuals.

The last two are described in Scenario Generation.

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

Externally-fitted PAR models

Novomodelo accepts fully user-supplied parameters — you may fit a PAR($p$) model with any method outside Novomodelo and hand it the two parquet files (the UserProvidedNoHistory and UserProvidedAll paths on the Configure tab). Novomodelo 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 Novomodelo’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$: Novomodelo 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, Novomodelo’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 Novomodelo 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: Novomodelo reconstructs $\psi_{m,\ell}$ and $\sigma_m$ from the stored std_m3s, so a mismatch is silent.
• Novomodelo 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 at or below a numerical floor 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; internally-fitted coefficients are not tested for stationarity but, like every load path, are guarded by the closure derivation’s own backstop — a non-finite derived $r_m$ is a hard load error naming the hydro and season.
• 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 Format and Configuration) — this applies on every load path, including a plain file-load with no fitting.
• The lags advance only under season_definitions. A per-stage season_id satisfies the rule above, but without season_definitions no stage completes a lag period and the derived inflow-lag seed is all zero, whatever inflow_history and recent_observations hold: the load succeeds (with a warning that the initial inflow lags seed to 0), and every stage receives those zero lags.

Outputs

The entire output/stochastic/ directory is written only when the case sets exports.stochastic: true (in the exports config block); it is off by default, so a default run produces no stochastic/ directory at all. The per-file conditions in the table below are secondary — they decide which artifacts appear inside the directory once stochastic export is enabled.

File	Role
output/stochastic/inflow_seasonal_stats.parquet	The seasonal statistics the run used, fitted or read from scenarios/, same schema as the input file — written on every stochastic export.
output/stochastic/inflow_ar_coefficients.parquet	The AR coefficients the run used, fitted or read from scenarios/, same schema as the input file — written on every stochastic export.
output/stochastic/inflow_annual_component.parquet	The PAR(p)-A annual triple the run used, fitted or read from scenarios/ — written on every stochastic export; it carries rows only for the stage models with an annual term (a classical PAR(p) model has none), so it can hold no rows.
output/stochastic/fitting_report.json	Per-hydro summary of the selected AR order and coefficients, a human-readable complement to the parquet outputs above — written only when Novomodelo estimated part or all of the inflow model (its seasonal statistics, its AR coefficients or both) from inflow_history.parquet.
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. To replay it on a subsequent run, follow Running Studies — Round-trip workflow, which covers the copy into scenarios/noise_openings.parquet and the openings declaration; the copy loads only when the case’s declared stage ids are 0, 1, 2, … (the export numbers stages by 0-based position, and the load refuses other ids) — 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 Novomodelo 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 Novomodelo always derives from the AR coefficients via the periodic-ACF closure (methodology §4). 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 (§4.1 of the methodology) — 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.

Residual-variance rejection (r_squared)

The closure sets $r_m^2 = 1 - \sum_\ell \psi^*_{m,\ell}\rho_m(\ell)$, the standardized residual (innovation) variance of season $m$ implied by $\psi^*$ (§4.1 of the methodology), and treats it as an admissibility test, not a quality score. Two guards act on it at load, the checks of methodology §5 invariant 2: each rejects the case with an error that names the hydro, and neither is a score on an accepted model. Braces in the quoted messages mark values the run fills in.

Supplied coefficients: the stationarity gate. A case that supplies inflow_ar_coefficients.parquet is checked by the PAR stationarity gate, which rejects a season whose $r_m^2$ is at or below 1e-12 (its other checks are listed in the Inputs & Outputs tab):

Hydro {hydro_id} season {season}: PAR stationarity gate rejected inflow_ar_coefficients.parquet -- implied residual variance r² = {r_squared} is at or below the numerical floor

Every load path: the derivation guard. After the model is assembled, the closure derivation of $r_m$ runs whatever the origin of the coefficients. It rejects the load when the closure system is singular, and when a derived $r_m$ is non-finite, which a negative implied variance $r_m^2$ produces:

residual_std_ratio closure is singular for hydro_id={hydro_id}

derived residual_std_ratio is non-finite for hydro_id={hydro_id} season={season}: the coefficients imply a non-stationary process (implied residual variance 1 - sum(psi*rho) is negative)

In the gate and non-finite messages season counts from 0: on a single-resolution season map it is the season’s position in calendar order; on a map that layers resolutions, or without a map, it is the rank of the season’s id among the season ids of the case. It equals the season id when the ids run $0, 1, 2, \ldots$ in calendar order.

A model fitted from the history is not gated and has no floor: only the derivation guard applies to it, so a fitted season whose $r_m^2$ is small but positive is accepted.

Novomodelo surfaces r_squared only inside the gate message: there is no separate CLI field or Python attribute for it, and no standalone variance signal on an accepted fit.

Order-selection behavior

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

• min_observations_per_season: a recommended minimum that only draws a warning, never an error; see Configuration — estimation for how groups are counted and when the warning is drawn.
• A season with a single observation rejects the case. Whenever Novomodelo estimates from inflow_history.parquet, it first estimates the seasonal statistics from the history, even when the case also supplies inflow_seasonal_stats.parquet (the supplied statistics then feed only the assembled model, not the fit). A (hydro, season) group with exactly one observation stops estimation with an insufficient-data error (”… need at least 2 for std estimation”), whatever min_observations_per_season is set to. A season with no observation has no statistics and is fitted at order 0.
• Automatic order reduction is logged per season in fitting_report.json (contribution_reductions[].reason), with three distinct reasons:
  • magnitude_bound — a prepass check: if any fitted coefficient exceeds max_coefficient_magnitude in absolute value, the season is reset to white noise (order 0, $r_m = 1$), independently of the contribution check.
  • phi1_negative — a prepass check, repeated after each re-fit: if the first fitted coefficient $\psi^*_{m,1}$ is negative, the season is likewise reset to white noise (order 0).
  • negative_contribution — the iterative Maceira–Damázio contribution check (methodology §3.6): the season’s order ceiling (initially max_order) drops by one, and PACF selection and the Yule-Walker fit re-run at the new ceiling, repeating until no season fails (a season whose ceiling reaches 0 takes order 0). The reduced_order recorded for the event is the largest order with no negative composed influence.

History classification runs before any order selection

Before PACF selection sees a (hydro, season) bucket, the history classification pipeline (methodology §3.3) 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 §7. 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 (§7.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 §7.7).

Inflow lags outside every cut

Whether a stage’s cuts carry the inflow-lag dimensions is the per-stage stages[].state_variables field. Its default, its validation warning, and whether that warning applies to supplied AR coefficients or to a model estimated from the inflow history are documented in Cut Management, the key’s single home.

Determinism guarantee

PAR(p) fitting is one of the model families covered by Novomodelo’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 & Provenance §3 for the full mechanism, which this chapter’s estimation pipeline shares with the computed-FPHA fitting pipeline described in Hydro Production Function Models.

Cross-References

• LP Formulation — The realized-inflow rows of the stage LP
• State Augmentation — AR inflow dynamics: state expansion, lag column pinning, reduced-cost extraction
• Inflow Non-Negativity Solution Methods — Methods for handling negative realizations produced by the PAR(p) model
• Scenario Generation — When external inflow scenarios are used in training, the backward pass turns the opening tree’s noise into inflows with the applied PAR model, in which hydros of AR order 0 and no annual component take the external samples’ mean and standard deviation.
• Notation Conventions — Defines inflow symbols ($a_{h,t}$, $\mu_m$, $\psi_{m,\ell}$, $\sigma_m$) and unit conventions