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.

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 {\epsilon }_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)}$: Residual standard deviation for season $m(t)$ (computed at runtime — see §3)
• ${\epsilon }_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,\dots ,M\}$ (e.g., $M=12$ for monthly, $M=52$ for weekly), the complete PAR(p) model requires:

3. Stored vs. Computed Quantities

The data model stores seasonal sample statistics and standardized AR coefficients with an explicit residual fraction. The relationship between stored and computed quantities is:

Stored in input files

These are provided in inflow_seasonal_stats.parquet and inflow_ar_coefficients.parquet (see Input Scenarios §3.1–3.2):

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 }^{\ast }$ is the direct output of the Yule-Walker fitting procedure (see §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 }^{\ast }\cdot \frac{s_m}{s_{m-\ell }}$$

Computed at runtime

From the stored quantities, 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 }^{\ast }\cdot \frac{s_m}{s_{m-\ell }}$$

Residual standard deviation (for noise scaling):

$${\sigma }_m=s_m\cdot {\text{residual\_std\_ratio}}_m$$

No autocorrelation values are needed at runtime. All required quantities are derived solely from the stored seasonal stats and AR coefficient file.

Why store residual_std_ratio rather than ${\sigma }_m$ directly? The residual std decomposes as ${\sigma }_m=s_m\cdot {\text{residual\_std\_ratio}}_m$, where $s_m$ is a conditioning quantity (swappable for climate scenario studies) and the ratio is a model dynamics property (fixed per PAR fit). Storing ${\sigma }_m$ directly would bake in a specific $s_m$: when the user swaps seasonal stats for a different climate scenario, the stored ${\sigma }_m$ would be stale and noise scaling would be inconsistent with the new variability level. Storing the ratio preserves correct proportionality — if seasonal variability changes, noise scales proportionally. See also PAR Coefficient Storage design document §3.4.

LP coefficients

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

$$a_h=\underset{\text{deterministic base}}{\underbrace{\left({\mu }_m-\sum _{\ell =1}^p{\psi }_{m,\ell }{\mu }_{m-\ell }\right)}}+\underset{\text{lag contribution}}{\underbrace{\sum _{\ell =1}^p{\psi }_{m,\ell }\cdot a_{h,\ell }}}+\underset{\text{stochastic innovation}}{\underbrace{{\sigma }_m\cdot {\eta }_t}}$$

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

4. Model Order Selection

The PAR order $p$ can vary by season. Available selection criteria:

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,\dots ,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 (§5.4) at order $k$. The last coefficient ${\hat{\psi }}_{m,k}^{\ast }$ 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 (§5.4).

Post-selection validation: After PACF selection, two rejection gates are applied iteratively:

• Negative ${\varphi }_1$ rejection: If ${\hat{\psi }}_{m,1}^{\ast }<0$ (first AR coefficient is negative), the order is reduced. Negative ${\varphi }_1$ contradicts the hydrological persistence property (inflows are positively autocorrelated at lag 1).
• Contribution-based validation: The recursively-composed contributions for each lag are computed. If any contribution is negative (indicating potential model instability), the order is reduced to the maximum lag with non-negative contributions. This implements NEWAVE’s reducao_ordem algorithm.

The reduction process is iterative: after each reduction, the PACF selection and coefficient estimation are re-run at the new ceiling, and the validation checks are repeated until all seasons pass or reach order 0.

4.2 AIC (Akaike Information Criterion)

$${\text{AIC}}_m(p)=N_m\ln ({\hat{\sigma }}_m^2)+2p$$

4.3 BIC (Bayesian Information Criterion)

$${\text{BIC}}_m(p)=N_m\ln ({\hat{\sigma }}_m^2)+p\ln (N_m)$$

4.4 Coefficient Significance

Include lag $\ell$ only if $|{\hat{\psi }}_{m,\ell }|>2/\sqrt{N_m}$.

In all methods, $N_m$ is the number of historical observations for season $m$.

5. Fitting Procedure

This section documents the five-step procedure for fitting PAR(p) parameters from historical inflow data. The fitting is performed when the system derives parameters from inflow_history.parquet (see Input Scenarios §2). When pre-computed parameters are provided directly in inflow_seasonal_stats.parquet and inflow_ar_coefficients.parquet, this procedure is not executed.

5.1 Notation

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

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-1}\sum _{t:m(t)=m}(a_{h,t}-{\bar{a}}_m{)}^2}$$

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-1}\sum _{t:m(t)=m}\left(a_{h,t}-{\bar{a}}_m\right)\left(a_{h,t-\ell }-{\bar{a}}_{m-\ell }\right)$$

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}^{\ast },\dots ,{\psi }_{m,p}^{\ast }$ 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: For row $i$ and column $j$ (0-indexed, $0\le i,j<p$), the reference season is shifted by row index:

$$[{\mathbf{R}}_m{]}_{i,j}={\hat{\rho }}_{(m-i)\mathrm{mod}M}(|j-i|)$$

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^{\prime }}(0)=1$ for any season $m^{\prime }$). The matrix is symmetric but not Toeplitz when $M>1$, because each row references a different season for its autocorrelation values.

RHS construction: Each RHS element also uses a shifted reference season:

$$[{\mathbf{r}}_m{]}_i={\hat{\rho }}_{(m-i)\mathrm{mod}M}(p-i)$$

This comes from column $p$ of the extended $(p+1)\times (p+1)$ version of the periodic autocorrelation matrix.

The full system is:

$$\left(\begin{array}{ccccc}1& {\hat{\rho }}_m(1)& {\hat{\rho }}_m(2)& \cdots & {\hat{\rho }}_m(p-1)\\ {\hat{\rho }}_{(m-1)}(1)& 1& {\hat{\rho }}_{(m-1)}(1)& \cdots & {\hat{\rho }}_{(m-1)}(p-2)\\ {\hat{\rho }}_{(m-2)}(2)& {\hat{\rho }}_{(m-2)}(1)& 1& \cdots & {\hat{\rho }}_{(m-2)}(p-3)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\hat{\rho }}_{(m-p+1)}(p-1)& {\hat{\rho }}_{(m-p+1)}(p-2)& {\hat{\rho }}_{(m-p+1)}(p-3)& \cdots & 1\end{array}\right)\left(\begin{array}{c}{\psi }_{m,1}^{\ast }\\ {\psi }_{m,2}^{\ast }\\ {\psi }_{m,3}^{\ast }\\ ⋮\\ {\psi }_{m,p}^{\ast }\end{array}\right)=\left(\begin{array}{c}{\hat{\rho }}_m(p)\\ {\hat{\rho }}_{(m-1)}(p-1)\\ {\hat{\rho }}_{(m-2)}(p-2)\\ ⋮\\ {\hat{\rho }}_{(m-p+1)}(1)\end{array}\right)$$

where all season indices are taken modulo $M$.

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

where:

• ${\mathbf{R}}_m$ is the $p\times p$ periodic correlation matrix (symmetric but not Toeplitz for $M>1$)
• ${\mathbf{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{\mathbf{\psi }}}_m^{\ast }={\mathbf{R}}_m^{-1}{\mathbf{r}}_m$$

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

5.5 Step 4 — Store Standardized Coefficients and Residual Fraction

The Yule-Walker solution ${\psi }_{m,\ell }^{\ast }$ 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.

Compute and store the residual std ratio:

$${\widehat{\text{residual\_std\_ratio}}}_m=\sqrt{1-{\mathbf{\psi }}_m^{\ast \top }{\mathbf{r}}_m}=\sqrt{1-\sum _{\ell =1}^p{\psi }_{m,\ell }^{\ast }\cdot {\hat{\rho }}_m(\ell )}$$

Both ${\psi }_{m,\ell }^{\ast }$ (one row per lag) and ${\widehat{\text{residual\_std\_ratio}}}_m$ (repeated across all lag rows of the same (hydro, stage) group) are written to inflow_ar_coefficients.parquet.

5.6 Step 5 — Residual Standard Deviation

The residual standard deviation for season $m$ is recovered at runtime from the stored ratio (see §3):

$${\hat{\sigma }}_m={\hat{s}}_m\cdot {\widehat{\text{residual\_std\_ratio}}}_m$$

For reference, the full expression in terms of fitting quantities is:

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

6. Validation Invariants

After fitting or loading pre-computed parameters, the following invariants must hold:

1. Positive residual variance: ${\sigma }_m^2>0$ for all seasons. If violated, the AR model explains all variance — likely overfitting.
2. Stationarity: Roots of $1-{\sum }_{\ell }{\psi }_{m,\ell }{z}^{\ell }=0$ lie outside the unit circle. Ensures the AR process is stable and does not diverge.
3. Correlation matrix positive definite: ${\mathbf{R}}_m$ is invertible. Required for Yule-Walker solution to exist. If violated, the historical record may be too short for the requested AR order.
4. No systematic bias: Residuals ${\epsilon }_t$ have mean near zero. Indicates the model captures the mean structure correctly.
5. AR order derivation: The number of coefficient rows per (hydro_id, stage_id) in inflow_ar_coefficients.parquet determines the AR order $p_m$. Lags must be contiguous: $\{1,2,\dots ,p_m\}$.
6. Residual std ratio consistency: The residual_std_ratio value must be identical across all lag rows sharing the same (hydro_id, stage_id) group, and must lie in $(0,1]$.

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, scaled by the seasonal standard deviation. In fully standardized form:

$$\frac{a_{h,t}-{\mu }_{m(t)}}{{\sigma }_{m(t)}}=\sum _{\ell =1}^p{\varphi }_{m(t),\ell }\frac{a_{h,t-\ell }-{\mu }_{m(t-\ell )}}{{\sigma }_{m(t-\ell )}}+{\epsilon }_t$$

where:

• ${\varphi }_{m(t),\ell }$: AR coefficients in fully standardized form (correlations between normalized deviations)
• ${\sigma }_{m(t)}$: residual standard deviation for season $m(t)$ (derived at runtime — see §3)
• ${\epsilon }_t\sim \mathcal{N}(0,1)$: innovation noise

The input files store ${\psi }_{m,\ell }^{\ast }$ (standardized by seasonal std $s_m$, not residual std ${\sigma }_m$). The next step converts these to original-unit ${\psi }_{m,\ell }$ for use in the LP.

7.2 Coefficient Conversion

The stored standardized coefficients ${\psi }_{m,\ell }^{\ast }$ 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 }^{\ast }\cdot \frac{s_m}{s_{m-\ell }}$$

The residual standard deviation is also derived at this preprocessing step:

$${\sigma }_m=s_m\cdot {\text{residual\_std\_ratio}}_m$$

These conversions are performed once at LP construction time. They require only the seasonal stats ($s_m$) and the stored model quantities (${\psi }_{m,\ell }^{\ast }$, ${\text{residual\_std\_ratio}}_m$) — no autocorrelation values, no historical data.

7.3 LP-Ready Form

Multiplying both sides of the canonical form (7.1) by ${\sigma }_{m(t)}$ and rearranging yields the LP-ready equation:

$$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 {\epsilon }_t$$

where ${\psi }_{m(t),\ell }$ and ${\sigma }_{m(t)}$ are derived from stored quantities as described in §7.2.

This decomposes the inflow into three additive components:

1. Lag contribution: $\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: ${\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 {\epsilon }_t$ — noise draw scaled by the seasonal residual 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 (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 {\epsilon }_t$$

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 §5), 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 (see LP Formulation §5a):

$$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 {\epsilon }_t$$

where:

• $b_{h,m(t)}$ is read from PrecomputedParLp.deterministic_base[stage][hydro]
• ${\sigma }_{m(t)}$ is read from PrecomputedParLp.sigma[stage][hydro]
• ${\epsilon }_t$ is the scenario noise draw for this (stage, hydro)

The ${\psi }_{m(t),\ell }$ coefficients from PrecomputedParLp.psi[stage][hydro][lag] 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 precomputation in PrecomputedParLp eliminates all redundant arithmetic from the hot path.

7.6 Summary of LP Components

Cross-References

• Input Scenarios §3.1–3.2 — Defines inflow_seasonal_stats.parquet (μ, s) and inflow_ar_coefficients.parquet (ψ* per lag, residual_std_ratio)
• LP Formulation §5 — AR inflow dynamics in the LP: state expansion, lag fixing constraints, dual variables
• Internal Structures §14 — PrecomputedParLp struct caching the three LP components derived in section 7
• Inflow Non-Negativity — Methods for handling negative realizations produced by the PAR(p) model
• Scenario Generation §4.2 — 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 (§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