Martin Mashalov / Research
Energy Forecast Lab / Technical record / September 2026

The mathematics of Energy Forecast Lab

Implemented methods, retained evidence and the limits of the claim.

On this page · 7 sections

This project studies whether simple updates to an operator's electricity-demand forecast improve daily-peak predictions, then builds the recording and accounting needed to test their economic value. Historical forecast diagnostics and prospective paper execution are separate experiments. No historical forecast score is presented as a trading return.

The implementations are retained under their original weatherpred Python namespace for source identity; the documented experiments here concern energy.

1. Define the target before fitting

Let Dd,hD_{d,h} be hourly average system demand for operating day dd, and let Fd,hF_{d,h} be its published operator forecast. The historical proxy target is

Yd=maxhHdDd,h.Y_d=\max_{h\in H_d}D_{d,h}.

HdH_d is the exact America/Chicago operating-day grid: 23, 24 or 25 physical hours, identified by UTC hour endings. Midnight belongs to the hour's operating date, not automatically to the date printed on its ending timestamp. Energy in MWh over a full hour is numerically equal to average MW over that hour.

Historical experiments use raw Demand and Demand forecast from the EIA balancing-authority workbook. Adjusted or imputed demand is not silently substituted. Demand/forecast range, imputation and UTC-identity flags qualify labels; a generic unrelated quality flag is not sufficient grounds to exclude an observation. The hourly dataset retains missing and flagged cells rather than deleting whole days.

This is different from the inspected contract's target: the highest unrounded TOTAL in the first complete daily ERCOT NP6-346-CD Actual System Load by Forecast Zone CSV. Later revisions do not redefine that contract value. The originally proposed NP6-345-CD Weather Zone product is not interchangeable. Contract rules determine the operating day; a date-like ticker suffix can refer to the following day.

The EIA workbook was received in September 2026 and contains revised history, not a complete archive of first publications. Its modification time is not the original availability time of each cell. Historical studies explicitly assume forecast availability at noon Chicago, demand availability one hour after hour end, and label availability at day+2 noon. EIA's definitions and availability notes explain its reporting system; those notes do not verify every historical receipt assumed by this experiment.

2. F001: correct the operator's peak error

Define the raw forecast peak Bd=maxhFd,hB_d=\max_h F_{d,h} and residual ed=YdBde_d=Y_d-B_d. Three fixed methods predict:

Method Peak prediction
Operator Y^d=Bd\widehat Y_d=B_d
Mean residual correction Y^d=Bd+etrain\widehat Y_d=B_d+\overline e_{\mathrm{train}}
Ridge residual correction Y^d=Bd+zdβ\widehat Y_d=B_d+z_d^\top\beta

The ridge features are the exact seven-calendar-day-lagged peak residual, annual seasonal sine/cosine, and today's operator peak. Missing the exact lag makes the case unavailable; the code does not search for a convenient earlier replacement.

Training means and standard deviations transform features, with an intercept appended. For standardized design matrix ZZ:

β=(ZZ+Λ)1Ze,Λ=diag(0,10,10,10,10).\beta=(Z^\top Z+\Lambda)^{-1}Z^\top e,\qquad \Lambda=\operatorname{diag}(0,10,10,10,10).

The implementation solves the linear system rather than explicitly calculating the inverse. The intercept is unpenalized; the penalty of 10 is fixed, not selected from evaluation scores. Twelve monthly fits use the preceding 365 calendar days, with labels strictly available before each month's UTC cutoff and at least 180 usable training days. The same eligible evaluation days are used for all three methods.

Code: energy_baselines.py, F001 runner. Results: F001 evidence.

3. F002: update the remaining hourly path

At noon, 14:00 and 15:00 Chicago, only observations whose hour ending plus one hour is no later than the decision are visible. The remaining unresolved hours include the hour ending at the decision itself. No current-day future actual values or retrospective quality flags enter the prediction view.

Let u1,u2,u3u_1,u_2,u_3 be actual minus forecast demand for the last three available hours, and uˉ=(u1+u2+u3)/3\bar u=(u_1+u_2+u_3)/3. The fixed controls use remaining-hour centers

μd,h=Fd,h+wuˉ,w{0,12,1},\mu_{d,h}=F_{d,h}+w\bar u,\qquad w\in\{0,\tfrac12,1\},

named operator, half update and full update. Ridge instead learns a separate correction for each remaining hour:

μd,h=Fd,h+zdβh.\mu_{d,h}=F_{d,h}+z_d^\top\beta_h.

Its six fixed features are recent mean residual, last-minus-first residual, observed-peak error relative to the corresponding forecast prefix, remaining forecast peak, and seasonal sine/cosine. A single matrix solve fits all hourly response columns using penalty 10 and an unpenalized intercept. Features are standardized only on training data.

For evaluation month mm, training covers a trailing 365-calendar-day window ending at m63m-63 days. A disjoint 60-day residual library covers [m63,m3)[m-63,m-3). Its labels must be available before mm; training labels must precede the earlier fit cutoff. Actual eligible support was 363 training and 60 library days per month/clock.

Whole-day residual scenarios

For library day jj, save the entire remaining-hour error vector

rj,h=Dj,hμj,h.r_{j,h}=D_{j,h}-\mu_{j,h}.

For a new decision, each donor vector produces one scenario:

Xj=max(Od,maxhmax(0,μd,h+rj,h)),X_j=\max\left(O_d,\max_h\max(0,\mu_{d,h}+r_{j,h})\right),

where OdO_d is the observed maximum so far. Whole vectors preserve dependence between hours. Independently sampling each hour would change the distribution of the daily maximum. There are 60 donor days, not thousands of independent simulations.

The historical OdO_d comes from revised observations. Treating it as a fixed floor is an oracle assumption: provisional real-time observations can change, and EIA is not the exact settlement source. The forward study therefore adds explicit sensitivity and entry-exclusion rules described below.

Historical exceedance probability is

p^(K)=160j1{Xj>K}.\widehat p(K)=\frac1{60}\sum_j\mathbf1\{X_j>K\}.

It is monotone across increasing thresholds, with probability increments of 1/601/60 before ties. Thresholds are the operator peak plus fixed offsets 5000,2500,0,2500,5000-5000,-2500,0,2500,5000 MW. They are synthetic forecasting thresholds, not a census of executable market contracts. The operator distribution already includes donor residual calibration and the observed maximum; it is a stronger baseline than an unadjusted operator point forecast.

Code: energy_intraday.py, F002 runner. Results: F002 evidence.

4. F003: tune dispersion on a separate window

F002 revealed undercoverage. F003 changes the spread of whole donor vectors while preserving their hourly means:

rj=rˉ+s(rjrˉ),s{1,1.25,1.5,2}.r'_j=\bar r+s(r_j-\bar r),\qquad s\in\{1,1.25,1.5,2\}.

Scale 1 preserves the original vectors exactly. One common scalar applies to every hour of a vector. Hourly means are preserved, but the mean or median of the nonlinear daily maximum can change.

For each month, clock and method, the scale with the lowest mean 90% interval score on a separate tuning block is selected; ties choose the smaller scale. All four choices' tuning scores are retained. Both half update and ridge receive this opportunity, alongside unchanged controls and the operator baseline.

Role Calendar window relative to month start mm Availability cutoff
Fit Trailing 365 days ending at m126m-126 m126m-126, 00:00 UTC
Residual library [m126,m66)[m-126,m-66) m63m-63, 00:00 UTC
Tune scale [m63,m3)[m-63,m-3) mm, 00:00 UTC
Evaluate The following calendar month Models and scale already frozen

Every label must satisfy the strict assumed release gate for its role. These windows are disjoint, not statistically independent weather samples. Unchanged F003 controls use the same older windows as scaled variants. Comparing F002's unchanged model directly with F003's scaled model would mix dispersion with training recency.

Code: energy_dispersion.py. Frozen protocol, historical proposal, F003 evidence.

5. Score the point, distribution and tails separately

For predicted median X~d\widetilde X_d and actual peak YdY_d, report MAE, RMSE and signed bias. MAE is the mean of X~dYd|\widetilde X_d-Y_d|; RMSE is the square root of mean squared error. F001 uses its point predictions directly.

For equally weighted peak scenarios, the continuous ranked probability score is

CRPS=1njXjY12n2j,kXjXk.\mathrm{CRPS}=\frac1n\sum_j|X_j-Y| -\frac1{2n^2}\sum_{j,k}|X_j-X_k|.

It scores the whole distribution in MW, with smaller values better. Brier score averages (p^(K)1{Y>K})2(\widehat p(K)-\mathbf1\{Y>K\})^2 over the five fixed thresholds. A good average Brier score does not certify probabilities on the subset later chosen for trading.

Intervals use empirical linear-interpolated quantiles. At nominal coverage 1α1-\alpha, L=Qα/2L=Q_{\alpha/2} and U=Q1α/2U=Q_{1-\alpha/2}. Report coverage 1{LYU}\mathbf1\{L\le Y\le U\}, width ULU-L, and separate lower/upper miss rates. F003's primary 90% interval score is

IS90=(UL)+20(LY)1{Y<L}+20(YU)1{Y>U}.\mathrm{IS}_{90}=(U-L)+20(L-Y)\mathbf1\{Y<L\} +20(Y-U)\mathbf1\{Y>U\}.

This penalizes both excessive width and missed outcomes. Reaching 90% coverage by widening indefinitely would not optimize the score. F003 predeclared a 5% score improvement as a practical diagnostic threshold, separate from whether an interval for the mean difference excludes zero; it was not a test that improvement exceeds 5% and not a profit hurdle.

Paired calendar bootstrap

Each study retains all 365 evaluation dates, including missing cases. For every contrast it forms daily paired score differences, then resamples circular blocks of seven calendar days. Missing dates keep their locations. Shared indices are used across comparisons; each replicate is the sum of available differences divided by its available-day count. There are 10,000 draws, without rerolls. A zero-support draw makes its interval unavailable.

Study Fixed comparison family Seed Two-sided percentile interval
F001 Two candidate-versus-operator MAE differences 20260907 97.5% per contrast
F002 Three pairs × CRPS/Brier × three clocks = 18 20260908 α=0.05/18\alpha=0.05/18 per contrast
F003 Four pairs × interval score/CRPS × three clocks = 24 20260909 α=0.05/24\alpha=0.05/24 per contrast

F002's pairs are full–operator, half–full and ridge–full. F003's pairs are scaled ridge–ridge, scaled half–half, scaled ridge–scaled half and scaled ridge–half. The intervals use quantiles α/2\alpha/2 and 1α/21-\alpha/2; the family adjustment is Bonferroni. They remain descriptive: dependence assumptions, revised history and reuse of 2025 as a development period prevent an untouched-holdout claim.

6. T002: a fixed prospective paper experiment

T002 froze three methods—operator, half update and ridge—at all three clocks for September 7, 2026. F004 preparation used revised history actually received before future decisions, with 359 qualified training days and 60 separate calibration days. This is permissible for future predictions; it does not mean those models existed at earlier historical dates. F003 was not substituted into the frozen pilot.

Live input selection freezes the first single complete daily DF response received by noon. Later revisions do not replace that forecast. For each permitted D hour, it takes the latest eligible actually received raw reported revision. Missing or imputed revisions invalidate an earlier value rather than allowing fallback. UTC identity, exact day grids, body hashes, source identities and receipt gates are retained. Publication, local network receipt and persistence are distinct facts; see the availability evidence.

Side scores and integer sizing

For 60 scenarios, the central score uses additive smoothing:

pc(K)=1+j1{Xj>K}62.p_c(K)=\frac{1+\sum_j\mathbf1\{X_j>K\}}{62}.

The fixed sensitivity scores are

pY(K)=1+j1{Xj500>K}62,pN(K)=1+j1{Xj+500K}62.p_Y(K)=\frac{1+\sum_j\mathbf1\{X_j-500>K\}}{62},\qquad p_N(K)=\frac{1+\sum_j\mathbf1\{X_j+500\le K\}}{62}.

The 500 MW shifts are assumptions, not measured error bounds or confidence intervals. The two stressed side scores need not sum to one. Thresholds at or below the provisional observed peak plus 500 MW are excluded, preventing an observation floor from forcing an entry.

The ask comes from the opposite-side bid: aY=1bNa_Y=1-b_N, aN=1bYa_N=1-b_Y. For ask aa, the original limit is =a+0.020.01\ell=\lceil a+0.02\rceil_{0.01}. Let MM be the verified fee multiplier and define the per-contract planning allowance

A=+0.0175M+0.02.A=\ell+0.0175M+0.02.

Require pA0.10p-A\ge0.10 and <1\ell<1. For cash CC, risk fraction ρ{0.05,0.20}\rho\in\{0.05,0.20\} and observed whole-contract depth QQ, enumerate integer quantities

0qmin(10000,Q,ρC/A)0\le q\le\min\left(10000,Q,\left\lfloor\rho C/A\right\rfloor\right)

and maximize

g(q)=plogCAq+qC+(1p)logCAqC.g(q)=p\log\frac{C-Aq+q}{C}+(1-p)\log\frac{C-Aq}{C}.

Cash, q=0q=0, is a genuine alternative. Candidate ties have deterministic ticker/side ordering. This is conditional log-wealth optimization using an experimental proxy score, not proof that pp is the true probability or that this is an optimal real trade.

Each account permits one entry attempt per daily event, including an unfilled attempt. All thresholds and clocks share its budget. Three methods × two risk fractions × three execution scenarios give 18 alternative $200 accounts, not $3,600 deployed or 18 independent outcomes.

Delayed conditional fills

An immutable intent precedes execution inspection. Arrival is 60 seconds after intent; expiry is 180 seconds after intent. The first relevant later recorded request determines the attempt. Its request, receipt and processing must fit the window, with at most ten seconds of response time. A bad first request cannot be replaced by a favorable later book.

Scenario Displayed depth retained Price stress, once Account precision
Optimistic 100% $0.00 $0.0001
Base 50% $0.01 $0.0001
Pessimistic 25% $0.02 $0.01

Walk opposing bids in price order; each implied ask plus stress must remain within the original limit. Retained depth is floored to 0.01 contracts at each level. Initial orders have integer quantities, while a conditional partial fill can be fractional. A displayed quote or touched price does not establish actual execution; passive queue fills are not inferred.

Fees, cash and unresolved exposure

For a taker slice of quantity qq at price pp, the supported quadratic model is

fmodel=0.07Mqp(1p).f_{\mathrm{model}}=0.07Mqp(1-p).

Round that trade fee upward to 10610^{-6} dollars. For account grid ε\varepsilon, let R=pqR=-pq for a buy, ff be the rounded trade fee, and cc the order's prior rounding carry:

b=Rfε,u=Rfb,v=min(c+uε,f+uε).b=\lfloor R-f\rfloor_{\varepsilon},\quad u=R-f-b,\quad v=\min(\lfloor c+u\rfloor_{\varepsilon},\lfloor f+u\rfloor_{\varepsilon}).

The cash change is b+vb+v, net fee is f+uvf+u-v, and new carry is c+uvc+u-v. A separate accumulator belongs to each order. This follows the inspected fee-rounding documentation; current series parameters and event overrides must also be supported. The shared fee helper supports maker coefficients, but T002 executes only taker scenarios.

The planning allowance is not a universal fee ceiling: many small fragments with cent balance rounding can exceed it. An overflow, missing first snapshot or unresolved execution locks the reserve and makes performance incomplete. Only a valid observed book can establish a conditional zero-fill. Unknown execution is not a free cancellation.

Cash is reserved at intent, debited for supported fills, and released from inventory only after an actual receipt proves a consistent finalized payout and any known later settlement time has arrived. Repeated final receipts cannot pay twice. Cancellation, refund, changed semantics or conflicting outcomes remain unresolved; they are not silently assigned an ordinary win/loss. EIA values do not determine paper payouts.

Code: energy_paper.py, trade inputs, fees, independent fill audit. See T002 protocol and audit protocol. The preserved protocols are historical frozen records, not instructions to extend their fixed collection windows.

7. Audit scope and interpretation

Independent routines reconstruct cash, exact fill arithmetic and integer selection from retained evidence. Forecast reconstruction in the operational auditor reuses frozen input/model functions, so it is not a second statistical model implementation. F003's separate auditor independently recalculated saved-coefficient paths, scaling, scores and all bootstrap contrasts without refitting models.

Source hashes and local append-only chains detect changed bytes and preserve ordering; they are not external notarization or proof of original publication. Synthetic tests verify implementation behavior, not strategy performance. A profitable-looking conditional forecast, a successful arithmetic audit, and complete economic evidence are three different claims.