PAR(p) Autoregressive Models

What Is a PAR(p) Model?

A Periodic Autoregressive model of order p (PAR(p)) is a time series model designed for data with strong seasonal patterns. It extends the classical autoregressive (AR) model by allowing every parameter to vary by season — the coefficients that govern January inflows are different from those that govern July inflows.

The “order p” indicates how many past time steps the model looks back. A PAR(3) model for a given month predicts the current inflow using the inflows from the previous three months. The order can differ by season: January might need only one lag while April might need four, reflecting different hydrological dynamics across the year.

The PAR(p) Equation

For hydro plant $h$ at stage $t$ falling in season $m(t)$, the PAR(p) model is:

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

In words: the inflow at stage $t$ equals the seasonal mean, plus a weighted combination of how much recent inflows deviated from their seasonal means, plus a random shock.

Parameters by Season

Each season $m$ has its own set of parameters:

The seasonal mean ${\mu }_m$ and sample standard deviation $s_m$ are estimated from historical data. The AR coefficients ${\psi }_{m,\ell }$ are fitted using the Yule-Walker equations (see below). The residual standard deviation ${\sigma }_m$ is derived at runtime from the other parameters (it is not stored independently).

How Lags Become State Variables

In the SDDP framework, decisions at each stage depend on a set of state variables that summarize everything the optimizer needs to know from the past. For the PAR(p) model, the state variables are the lagged inflows:

$$\text{State at stage }t:(v_{h,t},a_{h,t-1},a_{h,t-2},\dots ,a_{h,t-p_{\max }})$$

where $v_{h,t}$ is the reservoir volume and $a_{h,t-\ell }$ are the lagged inflows needed by the autoregressive equation. Each lag adds one state variable per hydro plant to the SDDP subproblem.

This is significant for problem size: a system with 150 hydro plants and a maximum PAR order of 6 adds up to $150\times 6=900$ state variables beyond the reservoir volumes. The LP formulation includes constraints that “shift” lagged inflows forward from one stage to the next, ensuring the autoregressive structure is respected across the Bellman recursion.

Stored vs. Computed Quantities

Cobre stores the natural outputs of the fitting process:

• Stored: seasonal means (${\mu }_m$), seasonal sample standard deviations ($s_m$), AR order ($p_m$), and AR coefficients in original units (${\psi }_{m,\ell }$)
• Computed at runtime: the residual standard deviation ${\sigma }_m$, derived from the stored quantities to guarantee consistency

This design avoids redundancy — ${\sigma }_m$ is fully determined by the other parameters and recomputing it is inexpensive.

Yule-Walker Fitting Procedure

When fitting PAR(p) parameters from historical inflow data, the AR coefficients are estimated by solving the Yule-Walker equations — a linear system that relates the autocorrelations of the data to the model coefficients. The procedure has five steps.

Step 1 — Seasonal Statistics

For each season $m$, compute the sample mean and standard deviation from historical observations $\{a_{h,t}:m(t)=m\}$:

$${\hat{\mu }}_m=\frac{1}{N_m}\sum _{t:m(t)=m}{a}_{h,t}$$

$${\hat{s}}_m=\sqrt{\frac{1}{N_m-1}\sum _{t:m(t)=m}{\left(a_{h,t}-{\hat{\mu }}_m\right)}^2}$$

where $N_m$ is the number of historical observations for season $m$.

Step 2 — Seasonal Autocorrelations

Compute the cross-seasonal autocorrelation at lag $\ell$ for season $m$. The cross-seasonal structure arises because lag $\ell$ at season $m$ reaches back to season $m-\ell$ (cyclically):

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

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

Note that ${\hat{s}}_{m-\ell }$ is the standard deviation of season $m-\ell$, not of season $m$. This is the defining feature of a periodic (as opposed to stationary) autoregressive model.

Step 3 — Yule-Walker System

For each season $m$, the coefficients in standardized form ${\psi }_{m,1}^{\ast },\dots ,{\psi }_{m,p}^{\ast }$ satisfy:

$${\mathbf{R}}_m{\mathbf{\psi }}_m^{\ast }={\mathbf{r}}_m$$

where:

$${\mathbf{R}}_m=\left(\begin{array}{cccc}1& {\hat{\rho }}_{m-1}(1)& \cdots & {\hat{\rho }}_{m-p+1}(p-1)\\ {\hat{\rho }}_m(1)& 1& \cdots & {\hat{\rho }}_{m-p+2}(p-2)\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{\rho }}_m(p-1)& {\hat{\rho }}_{m-1}(p-2)& \cdots & 1\end{array}\right),{\mathbf{r}}_m=\left(\begin{array}{c}{\hat{\rho }}_m(1)\\ {\hat{\rho }}_m(2)\\ ⋮\\ {\hat{\rho }}_m(p)\end{array}\right)$$

The solution is:

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

The matrix ${\mathbf{R}}_m$ is not a standard Toeplitz matrix (because consecutive rows use different seasons’ correlations), but it has a similar structure. The correlation matrix must be positive definite for the solution to exist; if not, the historical record may be too short for the requested order.

Step 4 — Residual Standard Deviation

After solving the Yule-Walker system, the residual standard deviation for season $m$ is:

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

This equals ${\hat{s}}_m$ times the square root of the unexplained variance fraction. If ${\hat{\sigma }}_m^2\le 0$, the model overfits — it explains all historical variance, leaving no room for the noise term.

Step 5 — Convert to Original Units

The Yule-Walker solution yields coefficients in standardized form ${\psi }_{m,\ell }^{\ast }$ (dimensionless, relating standardized deviations). The LP requires original-unit coefficients:

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

These are computed once at initialization and used directly as LP constraint matrix entries.

Key Properties

• Periodicity: All parameters vary by season, matching the strong seasonality of hydrological data.
• Parsimony: The model order $p$ is selected per season using AIC, BIC, or coefficient significance tests, avoiding unnecessary lags.
• Stationarity: Fitted models are validated to ensure the AR process does not diverge — the characteristic polynomial roots must lie outside the unit circle.
• Positive residual variance: After fitting, ${\sigma }_m^2>0$ must hold for all seasons. A zero or negative residual variance indicates overfitting.

Further Reading

• Stochastic Modeling — overview of why inflow uncertainty matters and how PAR(p) fits into the SDDP workflow
• Inflow Non-Negativity — how negative PAR(p) realizations are handled
• PAR Inflow Model Specification — complete mathematical specification with unit conversion formulas, runtime preprocessing steps, and all validation invariants

See also: For implementation details and usage, see the Stochastic Modeling guide in the software book.