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 six steps.

Implementation status: As of v0.1.1, this full fitting procedure is implemented in cobre-stochastic’s estimation module. Steps 1–5 are carried out by the seasonal statistics, autocorrelation, and AR coefficient estimators; Step 6 selects the model order via AIC before the final coefficients are computed.

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.

Levinson-Durbin recursion: In the implementation, the Yule-Walker system is solved using the Levinson-Durbin recursion rather than direct matrix inversion. This algorithm exploits the near-Toeplitz structure of ${\mathbf{R}}_m$ to solve the system in $O(p^2)$ time and $O(p)$ space, compared to $O(p^3)$ for a general LU factorization. For the per-season orders typical in hydro studies ($p\le 12$), the difference is small in absolute terms, but the recursion also produces intermediate partial autocorrelation coefficients that are useful for diagnosing model fit.

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.

Step 6 — Model Order Selection (AIC)

Before Steps 3–5 are applied at the final model order, the implementation selects the order $p_m$ for each season $m$ using the Akaike Information Criterion (AIC). For a candidate order $p$, the AIC is:

$$AI{C}_p=N\ln \left({\hat{\sigma }}_p^2\right)+2p$$

where $N$ is the number of historical observations for the season and ${\hat{\sigma }}_p^2$ is the residual variance obtained by solving the Yule-Walker system at order $p$. The first term penalizes poor fit (large residual variance); the second term penalizes model complexity (many lags). The selected order minimizes the criterion over the candidate range:

$$p_m^{\ast }=\arg \underset{p\in \{0,1,\dots ,p_{\max }\}}{\min }AI{C}_p$$

where $p_{\max }$ is the user-specified maximum order. Order $p=0$ corresponds to a white-noise model (no autoregressive structure); its residual variance equals the sample variance ${\hat{s}}_m^2$.

Implementation: This criterion is implemented in select_ar_order_aic in the cobre-stochastic estimation module. The function evaluates $AI{C}_p$ for each candidate order and returns the minimizer. BIC and coefficient significance tests are recognized alternatives but are not implemented in v0.1.1.

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 the AIC criterion (implemented in v0.1.1 via select_ar_order_aic). BIC and coefficient significance tests are recognized alternatives but are deferred.
• 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.