Forecast Evaluation and Combinations

Module 7 · Scoring, testing, and averaging forecasts

Gary Cornwall

Econ 6376 · The George Washington University

Module 6 recap: a forecast is a point, an interval and a density

Point

\tilde y_{t+h\mid t}=\mathbb{E}\!\left[y_{t+h}\mid\mathcal{F}_t\right]

the headline number

Interval

\tilde y_{t+h\mid t}\pm c_h

with a stated coverage

Density

f(y_{t+h}\mid\mathcal{F}_t)

the whole distribution; the fan chart draws it

Uncertainty grows with the horizon.

$\text{Var(}e_{t+h\mid t}\text{)}=\sigma^2\sum_{i=0}^{h-1}\psi_i^2$

h=1

\sigma^2

. Stationary:

\to\gamma_0

. Unit root: keeps growing.

[About 2 minutes. Recap 1 of 4; everything on screen at once, no clicks.] Notation, say it once: Module 6 wrote ŷ_{t+h|t}; from here on forecasts get a tilde, hats mean fitted/estimated. Three objects from Module 6. The point forecast is the conditional mean, and that is a consequence of squared-error loss, not a modelling choice (next slide). The interval is the point plus or minus a multiple of the forecast-error standard deviation, with a stated coverage. The density forecast is the whole conditional distribution; the fan chart draws it (Module 6’s AR(1), phi = 0.7, sigma^2 = 1: the fan opens and stops at the unconditional band, dashed). Why the fan widens, in words: each extra step adds a shock not yet seen, so the variance is sigma^2 times the sum of squared psi-weights. sigma^2 at one step (the innovation variance we have carried since Module 1), converging to gamma_0 for a stationary model; for the d = 1 SARIMA on UNRATE the sum never converges. The psi-weighted error sum itself stays in the Module 6 notes; do not re-derive it here. Bridge to today: every forecast in Module 7 is one step ahead, h = 1, so the error variance we care about is the one-step one, sigma^2 for a correct model. [Sources] Module 6 notes 6.2 (point, interval, density; the conditional mean under MSE), 6.5 (psi-weight form, forecast-error variance, AR(1) closed form, fan chart). Figure: images/ar1_fan.png, copied from slides/module_06/images/ (made by slides/module_06/make_figures.R, section 5).

Module 6 recap: the loss picks the target

Under MSE, the best point forecast is theconditional mean

Under MAEconditional median

Under asymmetric linear loss with weight

\alpha

the

\alpha

-quantile

Huber and MAPE are losses too; this course scores with them rather than deriving a target. MAPE is scale-free but undefined at

y_t=0

Choose the loss before you see the results; compare within a loss, never between: same loss, same

y

, same test observations

N

Today’s DM test takes a loss as an input. Which one, and who chooses?

[About 2 minutes. Recap 2 of 4.] The loss function encodes what it costs to be wrong, and the cost picks the shape: explosive per-unit cost is MSE, constant per-unit cost is MAE, Huber is the compromise, a directional cost is the asymmetric loss (the Texas grid). The optimal forecast follows the loss: the conditional mean under MSE, the median under MAE, the alpha-quantile under asymmetric linear loss. Huber and MAPE are losses we score with in this course; we do not pair them with an optimal forecast. MAPE: scale-free but undefined at zero and asymmetric by construction; today’s MASE (section 7.3) is the fix. g, not L, for the loss: L is the lag operator. The prompt sets up section 7.4: Diebold-Mariano compares the loss of two forecasts, so the loss is chosen first, by you, and dm.test()’s power argument encodes it (power = 2 squared, power = 1 absolute). [Sources] Module 6 notes 6.7 (“The Concept”, “How the Cost Structure Maps to the Loss Function”, “The Formulas”: MSE -> conditional mean, MAE -> conditional median, asymmetric -> alpha-quantile; MAPE caveats; within, never between). Layout after Module 6’s “The loss picks the target” slide (slides/module_06/parts/_loss_functions.qmd). Figure: images/loss_shapes.png, made by figures/make_loss_shapes.R (the curve logic of slides/module_06/parts/make_loss_figures.R, delta = 2, alpha = 0.75; the axis carries today’s tilde).

Module 6 recap: train first, then test, in time order

Train: estimate here, onceUNRATE 1960-01 … 2017-12

Test24 months

1960origin: December 2017 → 2019

Never a random split. Shuffling leaks the future into training: Biff with the sports almanac.

Module 6: one origin, many horizons. One fit, one path,

h=1,\ldots,24

Module 7: many origins, one horizon. Fixed coefficients, rolling data.

Today: many origins, one horizon. Same parameters, the information set rolls.

[About 2 minutes. Recap 3 of 4.] The split: everything before the cutoff trains, everything after tests. Random splits are wrong for time series because they put December in the training set and June in the test set; the Back to the Future analogy, Biff and the almanac. Left panel, Module 6’s design: stand at December 2017 and forecast 24 months in one shot. One fit, one path (orange); the band widens with the horizon while the actuals (grey) glide down. Its errors at h = 1 to 24 share the early shocks; keep that in mind, because section 7.4 explains why the DM test cannot be run on that path. Right panel, today’s design: the same SARIMA with the same coefficients, but each month forecasts only the next one. Each orange step starts at the month’s actual (the origin) and ends at the one-step forecast; when the next actual arrives it joins the information set. The parameters never move. The y-axes differ on purpose: on the left panel’s scale the 24-step band would flatten the staircase. Correction to say aloud: Module 6’s materials promised that Module 7 would re-estimate as the origin rolls. It does not, in Core. Module 7 estimates once and rolls only the information set; the two re-estimating schemes (expanding and fixed-width rolling windows) are named on “Three ways to run an evaluation” and run only in an optional Deeper Dive (tsCV()). [Sources] Module 6 notes 6.3 (“Why Random Splits Are Wrong”, “Temporal Splits”, “Two Evaluation Designs”), 6.8 Step 1 (ts_split(unrate, test_periods = 24)); Module 7 notes 7.1 (“The PS3 Design: Many Origins, One Horizon”). Figures: images/split_one_origin.png and images/split_many_origins.png, made by figures/make_split_panels.R from cached data/UNRATE.csv with the Module 5/6 model (ARIMA(1,1,2)(1,0,1)[12] fit once on 1960-2017): forecast(h = 24) on the left, fitted(Arima(full, model = fit)) on the right.

Why Choose Only One?

Open Release Date Radar

[About 15 minutes. The live wisdom-of-crowds exercise. Click the button; it opens the app in a new tab. Not embedded, on purpose: it is a live, networked page. The slide carries only the link; say the procedure aloud.] Say aloud, for each movie: (1) read the clues and WRITE DOWN your year, committing before anyone speaks; (2) call out your guess, I type each one in, and the histogram and crowd mean build live; (3) reveal: the actual date and the crowd’s error; note YOUR error too. Ask before the first reveal: will the crowd mean beat most of you, or only some of you? What the app is: a single page that runs in the instructor’s browser. There is no student login and no backend; guesses reach it only through the one year box on the left, which accepts any four-digit year from 1900 to 2030 and records it the moment the fourth digit is typed. So the room calls guesses out (or hands up written cards) and you type them. It loads a fresh session of five movies from TMDb: Back to the Future is always movie 1 (the course’s anchor film), then one random pick from a 1970-1990 pool, one from 1990-2010, two from 2011 on. Per movie: the poster and clues (genre, director, lead actors, tagline, runtime; any year in a clue is scrubbed). As guesses arrive the centre panel draws a histogram of years with a dashed “Mean” line; the right panel shows n, the mean year, the sample variance and the mean converted to an estimated date. Once there are five or more guesses a “95% Prediction Interval” row appears: mean plus or minus t times s times sqrt(1 + 1/n). Read it correctly if anyone asks: it is an interval for ONE MORE GUESS from this crowd, not for the release date. The spread of opinion is not the accuracy of the average. Reveal draws a green “Actual” line, prints the actual date, the crowd error in years, months and days, and a verdict badge. Next moves on; after movie 5 the button reads Summary: a table of each movie’s actual year, crowd mean, error and n, the average error, and “Best Guess / Worst Guess” callouts. Those callouts are the MOVIES where the crowd was closest and farthest, not individual students; the app does not track who guessed what. Reset clears the current movie’s guesses. Running it: insist on independence (write first, then call out); the exercise only works if the guesses are not herded. Back to the Future first is deliberate: the room knows it and the clues are easy, so it calibrates the procedure. Needs internet: the movie data and posters come from TMDb, and the chart library loads from a CDN. No-wifi fallback: do it on paper with five movies you name aloud; collect the guesses on the board and compute the mean by hand. The debrief does not depend on the app. [Sources] The app’s source, grycrnwll.github.io/tools/release-date-radar.html (instructor’s site repo, tools/release-date-radar.html), read 2026-10-12; planning/module_07_plan.md, “Analogies in play” (movie-release-year guesses, the class average beats most individual guesses).

Debrief: the crowd average is a forecast combination

Your guess of the release date is a forecast based on a mental model.
Each of your models has partially overlapping information sets: you can exploit this!
Each of you made a forecast $\tilde y^{(m)}$ of the release year, $m=1,\ldots,M$ . The crowd mean is

$\tilde y^{(c)}=\frac{1}{M}\sum_{m=1}^{M}\tilde y^{(m)}$

the simple average: every weight

w_m=1/M

, nothing estimated.

The SPF’s published forecast is the median across ~37 professional forecasters: forecast combination as standard practice.

[About 4 minutes.] The connection to section 7.5 in three sentences. A guess is a forecast; the crowd mean is the simple average, the zero-estimation-error benchmark of section 7.6, Method 1. It tends to beat most individuals because imperfectly correlated errors diversify, the portfolio argument made formal later today (two-forecast algebra). The figure: professionals do the same thing. Every grey dot is one forecaster’s call, made four quarters ahead, of the unemployment rate, plotted at the quarter it was aiming at; the blue line is the SPF’s published number, the median of those dots. The vertical spread at any quarter is disagreement: different models, different information sets, the same thing as the spread of guesses in this room. The Philadelphia Fed does not publish any one forecaster; it publishes a combination (a median rather than a mean, so one wild call cannot drag it). Then point at 2020. The orange spike towers over every dot aimed at 2020Q2: nobody in the panel saw the pandemic coming, so the median could not either. Diversification removes idiosyncratic error, the part where forecasters disagree; it cannot remove common error, the part they all share. That is “correlation eats the gain” (section 7.5), and the combo2 widget shows it. (The dots aimed at 2021 sit high because those surveys were taken in mid-2020, after the shock.) Optional spoken seed for section 7.8: a combination is only as good as its pool. On UNRATE one of our six tracks (the expanding mean of a series with a unit root) is broken, and the simple average gives it a sixth of the vote. Caution: do not claim the combination always beats every individual; it does not, and the notes never say so. [Sources] notes 7.5 (Concept: From Selection to Combination; the portfolio analogy; “Correlation eats the gain”), 7.6 Method 1 (simple average, w_m = 1/M, no hat), 7.7 (the index fund and its broken-member caveat), 7.8 (“The Honest Reading”); notation dictionary, Forecasting Notation (w_m, hat convention, y-tilde^(c)). Figure: images/spf_unemp.png, made by figures/make_spf_figure.R from data/spf/ (Philadelphia Fed SPF, UNEMP6 = the quarterly-average unemployment rate four quarters after the survey quarter; surveys 1999Q1 to 2026Q3; median 37 respondents per survey) and data/UNRATE_latest.csv (FRED, current vintage, quarterly averages).

Where we are: nothing new on the mixing board

$y_t=\alpha+\delta t+\sum_{j=1}^{p}\phi_j y_{t-j}+\sum_{l=1}^{q}\theta_l\epsilon_{t-l}+\epsilon_t$

Every slider is on, seasonal block included. Module 6 projected one fitted equation forward; Module 7 puts several on the table and referees between them.

Term	L4	L5	L6	L7
$\alpha$ , $\delta t$	On	On	On	On
$\phi_j y_{t-j}$ (all $p$ )	On	On	On	On
$\theta_l\epsilon_{t-l}$ (all $q$ )	On	On	On	On
Seasonal block $(P,D,Q)_s$	Off	On	On	On
What we do with it	Select	Diagnose	Forecast	Evaluate and combine $M$ forecasts

Today’s arc, three questions:

1. How do I score several forecasts?

U_2

, MASE

2. Is the difference real?Diebold-Mariano

3. Should I even pick one?combinations

[About 2 minutes. Part A, Evaluation, starts here (where-we-are to dm-unrate, about 43 minutes). The required “Where We Are” slide (design principles section 3), placed after the debrief because the opening is fixed.] Nothing new switches on today. Every slider is on, the seasonal block included, the same board as Modules 5 and 6. What changes is what we do with the fitted equation: Module 4 selected it, Module 5 diagnosed it, Module 6 used one fitted equation to project forward, and Module 7 puts several fitted equations, plus a few naive rules, on the table at once and referees between them. The arc in three questions, and the order of the evening. (1) How do I score several forecasts? Two relative metrics, Theil’s U_2 and MASE. (2) Is the difference between two forecasts real? The Diebold-Mariano test. That answers the question Module 6 left open. (3) Should I even pick one? Forecast combinations, and the surprising reason the simplest combination is so hard to beat; the radar exercise was the preview. Timing plan for the session (minutes per slide are in each slide’s notes): opening about 25 (recaps 6, radar 15, debrief 4); Part A about 43; Part B, Combinations, about 67; close 3. Total about 138, under the 150 budget. [Sources] notes 7.1 “The Module 6 Handoff” (the three-part arc) and “Where We Are on the Mixing Board” (the board table, without its always-on epsilon_t row); course_overview_15kft.md L7 (story, four-step arc).

The PS3 design: many origins, one horizon

Train: estimate every model here, once

y_1,\ldots,y_T

T=N_{\text{total}}-24

Weights window

T{+}1,\ldots,T{+}12

Holdout

T{+}13,\ldots,T{+}24

UNRATE 1960-01train to 2017-12 · weights 2018 · holdout 2019

Parameters are estimated once, on

y_1,\ldots,y_T

never again

The information set rolls.

\tilde y_{t\mid t-1}

uses actuals up to

y_{t-1}

h=1

throughout

Weights window,

T_w=12

dateswhere combination weights are estimated

Holdout,

N=12

dateswhere everything is scored

One fit per model; at each date, the training coefficients meet the newest data:

T_origin <- length(unrate) - 24  # Dec 2017
train_un <- upto(unrate, T_origin)
fit1 <- Arima(train_un, order = c(1, 1, 2),
  seasonal = list(order = c(1, 0, 1),
                  period = 12))
# one-step forecasts, coefficients frozen:
fitted(Arima(upto(unrate, T_origin + 24),
             model = fit1))

Subset a ts with window() (here upto()), never unrate[1:T_origin]: square brackets drop the monthly frequency.

[About 4 minutes. Built in full: the combination-puzzle widget’s T_w bubble links here.] The one design the whole module runs on; the midterm specified it and Problem Set 3 uses it. Restate it, do not assume it was built. Three blocks. Train, 1 to T with T = N_total - 24 (UNRATE 1960 to December 2017): every model is estimated here, once. The weights window, T+1 to T+12 (2018), T_w = 12 dates: Part B estimates combination weights here and nowhere else. The holdout, T+13 to T+24 (2019), N = 12 dates: everything, single tracks and combinations, is scored here. Two rules. Parameters are estimated once and never again. The information set rolls: the forecast of y_t uses actuals through y_{t-1}; when y_t arrives it becomes information for the forecast of y_{t+1}, but the coefficients do not move. Every forecast is one step ahead, h = 1. The idiom: Arima(series, model = fit) re-applies the training coefficients to a longer series without re-estimating, and fitted() returns at each date the one-step forecast built from information up to the previous date (Module 6 section 6.7). Contrast with Module 6 (previous recap): one origin, many horizons there; many origins, one horizon here. Different error series, different statistical properties (section 7.4). The window() trap, one sentence: y[1:T] silently drops the frequency; with model = fit the stored seasonal period survives, but any fresh fit on such a subset (auto.arima(), a seasonal Arima()) loses the season. The code block is illustrative; the canonical upto() lives in helpers/track_utils.R. [Sources] notes 7.1 “The PS3 Design: Many Origins, One Horizon”, “Loading UNRATE and Building the Tracks” (upto(), window() vs square brackets, Arima(…, model = fit) + fitted()); notation dictionary (T, N, T_w, N_total, h = 1 throughout Lecture 7).

Six forecast tracks

Track	On UNRATE	In the pool as
`arima_1`	ARIMA(1,1,2)(1,0,1) $_{12}$	the Module 5/6 model
`arima_2`	ARIMA(1,1,0)	a simpler alternative
`arima_3`	`auto.arima()` on train	the automatic pick
`locf`	$k=1$	the naive benchmark
`avg_k`	$k=3$	smooths sampling noise
`meanf`	$k=t-1$	the long-run mean

$\tilde y_{t\mid t-1}=\frac{1}{k}\sum_{j=1}^{k}y_{t-j}$

$k$ -period average forecasts, never “moving averages”: in this course that is the MA(

q

) block.

Error:

e_{m,t}=y_t-\tilde y^{(m)}_{t\mid t-1}

. Out of sample and observable: not

\epsilon_t

, not a residual.

UNRATE 2018-2019, actuals in black. meanf is off the chart, near 6.0%.

Why

k=3

for UNRATE, and not

k=24

[About 6 minutes. LIVE CODING beat: load the cached UNRATE, upto(), fit the three ARIMAs once, tracks_on(), then the alignment check. The figure is the result, already on the slide.] M = 6 forecasters, indexed m. Three ARIMA tracks: arima_1 is the Module 5/6 SARIMA; arima_2 a simpler credible alternative; arima_3 whatever auto.arima() picks on train (on this training sample ARIMA(1,1,3)(2,0,2)[12], richer than the Module 5 choice). Three k-period average forecasts: the mean of the previous k actuals. k = 1 is locf, last observation carried forward, the naive random walk; k = 3 is avg_k; k = t - 1 is meanf, the expanding mean of everything so far. Say it once and firmly: these are not “moving averages”. Moving average means the MA(q) block, theta_l epsilon_{t-l}, a different object. The error: y_t minus forecaster m’s one-step forecast. Out of sample, observable. Not the innovation epsilon_t, not an in-sample residual e_t. Every forecast carries a tilde from today (the fitted() values on a refiltered model are forecasts, so they get the tilde too). Fallback: have unrate, fits_un, F_w and F_h pre-run in the session before class, in case the live coding runs long; the figure already shows the result. Live code, in this order: unrate <- load_cached_fred(…) and window() to 1960-2019 (N_total = 720, T = 696, December 2017); fit1, fit2, fit3 on train_un; tracks_on(unrate, fits_un, w_dates) and the same for h_dates. The track builder re-applies the training coefficients with Arima(series, model = fit) and takes fitted(); kavg() averages the previous k actuals. Then the alignment check, before any metric: all.equal(F_h[, “locf”], as.numeric(unrate)[h_dates - 1]) returns TRUE. If locf on date t is not y_{t-1}, every number downstream is wrong. The picture: five tracks hug the actuals and differ by hundredths. meanf sits near 6.0% while unemployment runs below 4%: on an I(1) series the expanding mean is a sixty-year average of a level that has wandered. Plant it; it decides section 7.8. The prompt (Module 2 callback): UNRATE behaves like an I(1) series, so almost all the information about next month sits in the latest value. Averaging 24 months drags the forecast back to levels the series left long ago; three months only smooths the household survey’s month-to-month sampling noise. The choice of k is a statement about persistence. [Sources] notes 7.1 “Six Forecast Tracks” (the track table, the k-period-average definition, the error definition, why k = 3), “Loading UNRATE and Building the Tracks” (chunks load, upto, fit-three, track-builder, alignment, plot-unrate-tracks). Canonical upto(), kavg(), tracks_on() in helpers/track_utils.R. Figure: images/unrate_tracks.png, made by figures/make_figures.R from the notes’ plot-unrate-tracks data (tr_df).

The course simulation: one track is genuinely better

$(1-0.7L)(1-0.5L^{12})\,y_t=\alpha+\epsilon_t,\qquad \epsilon_t\sim N(0,1),\qquad \alpha=5(1-0.7)(1-0.5)=0.75,\qquad \mu=5$

Track	In the simulation
`arima_1`	ARIMA(1,0,0)(1,0,0) $_{12}$ , mean	true specification
`arima_2`	ARIMA(1,0,0), mean	misses the season
`arima_3`	ARIMA(0,1,1)	a needless difference
`locf`, `avg_k`, `meanf`	$k=1$ , $3$ , $t-1$	as on UNRATE

Each ARIMA estimated once, on 240 training months; then 12 weights-window and 12 holdout dates.

Built so the optimal weights are far from

1/M

: a method that finds arima_1 should beat the simple average.

The course realization (seed 6376): the end of training, then the weights window and the holdout.

[About 3 minutes. Built in full: the combination-puzzle widget’s “correctly specified model” bubble links here.] Every topic today also runs on one simulated DGP, so we can see what the tools do when we know the truth. A stationary seasonal AR: (1 - 0.7L)(1 - 0.5L^12) y_t = alpha + epsilon_t, with alpha = 5(1 - 0.7)(1 - 0.5) = 0.75, so the process mean is 5, and standard normal innovations. Lag-polynomial signs as always: minus signs in the AR polynomials; both roots outside the unit circle, so it is stationary. The six simulated tracks mirror the UNRATE pool. arima_1 is ARIMA(1,0,0)(1,0,0)_12 with a mean: the true specification, estimated. arima_2 is an AR(1) with a mean: it misses the seasonal block. arima_3 is ARIMA(0,1,1): a needless difference, a plausible student error. Then locf, avg_k with k = 3, and meanf. In the simulation arima_3 is a fixed specification rather than auto.arima(), so every Monte Carlo replication compares the same six forecasters. Each ARIMA is estimated once on 240 training months. Why it was built this way (say it now, it pays off in Part B): if equal weights were already optimal, every estimated weight would be pure noise and the simple average would win trivially. So one track is genuinely better and the optimal weights are far from 1/M. A weighting scheme that finds arima_1 should beat the simple average; the question for section 7.7 is whether it can find it with twelve observations. In the combination-puzzle widget, the faint dotted line is arima_1 alone: what finding the dominant track is worth. The figure: the seed-6376 realization, the last four years of training, then the weights window and the holdout. Seed 6376 is the course realization used in sections 7.3, 7.4 and 7.6 and Stage 1. [Sources] notes 7.1 “The Course Simulation” (the DGP, the six tracks, chunk course-dgp: course_dgp(), fit_sim(), set.seed(6376), n_train = 240); 7.7 Concept (“The DGP was built so the experiment is fair”). Figure: images/course_dgp.png, made by figures/make_figures.R from the notes’ y_sim.

Three ways to run an evaluation

Scheme	Parameters	Estimation window at origin $t$	Fits per model
Fixed-parameter one-step (PS3)	estimated once	$1{:}T$ ; information rolls to $t-1$	1
Expanding window	re-estimated at every origin	$1{:}(t-1)$	one per origin
Fixed-width rolling window	re-estimated at every origin	$(t-W){:}(t-1)$	one per origin

Module 7 builds only the first row. The other two re-estimate at every origin.

Re-estimation buys robustness to change at the cost of noisier parameters. Fixed-width windows forget old regimes on purpose.

Which scheme would you use if you suspected a break in 2008?

[About 2 minutes. Minimal slide; point to the notes.] The question: as the origin moves forward, do you re-estimate, and on which data? Module 7 builds only the first row: parameters estimated once on 1 to T, the information set rolling to t - 1, one fit per model, so differences between tracks are differences between models, not between estimation samples. The other two rows re-estimate at every origin: expanding (all data to t - 1) and fixed-width rolling (the last W observations). Module 6 named these two alongside its own fixed-origin design (the left panel of recap 3, which has a single origin and so is not a row here). The correction of Module 6’s re-estimation promise was said on that recap slide; here just point back to it: we use the first, and tsCV() runs the other two in an optional Deeper Dive that stays off the deck. The dot figure, one row per forecast origin, one dot per date. Navy is the estimation sample, light blue is data used as information only (not re-estimated), orange is the date being forecast. In the PS3 panel the navy block never moves and the light blue grows; in the other two the navy block itself grows (expanding) or slides (fixed width). The prompt: if you suspected a break in 2008, the fixed-width window is the scheme that lets pre-2008 data drop out of the estimation sample; the expanding window keeps it forever, and the PS3 design never re-estimates at all. The two sentences to remember. Re-estimation buys robustness to change at the cost of more computation and noisier parameters. Fixed-width windows forget old regimes on purpose: useful if the world changed, costly if it did not, because each fit uses fewer observations. Four-step deviation by design: a definitional bridge, no Math or Simulate step. [Sources] notes 7.2 (the scheme table, the plot-schemes figure, the two sentences, “It is the only scheme this module builds”); Module 6 notes 6.3 “Two Evaluation Designs”. Figure: images/eval_schemes.png, made by figures/make_figures.R from the notes’ plot-schemes data (dots), one row of panels instead of one column.

$U_2$ and MASE answer different questions

RMSE and MAE are in units of

y

; MAPE is a percentage. Is 0.1 points good? Two relative scores, two different questions.

$U_2$ : a race. Did model

m

beat locf on these holdout dates?

$U_2^{(m)}=\frac{\text{RMSE}^{(m)}}{\text{RMSE}^{(\texttt{locf})}}$

Both RMSEs over the same holdout dates,

t=T{+}13,\ldots,T{+}24

U_2<1

beat locf here

U_2>1

lost to it here

MASE (mean absolute scaled error): a yardstick. Model

m

’s holdout miss, relative to the naive’s typical one-step miss in training.

$\text{MASE}^{(m)}=\frac{\text{MAE}^{(m)}_{\text{holdout}}}{\text{MAE}^{(\texttt{locf})}_{\text{train}}}$

$\text{MAE}^{(m)}_{\text{holdout}}=\frac{1}{N}\sum_{t=T+13}^{T+24}\bigl\lvert y_t-\tilde y^{(m)}_{t\mid t-1}\bigr\rvert$

$\text{MAE}^{(\texttt{locf})}_{\text{train}}=\frac{1}{T-1}\sum_{t=2}^{T}\bigl\lvert y_t-\tilde y^{(\texttt{locf})}_{t\mid t-1}\bigr\rvert$

with

\tilde y^{(\texttt{locf})}_{t\mid t-1}=y_{t-1}

. Seasonal variant:

y_{t-s}

replaces

y_{t-1}

, over

t=s{+}1,\ldots,T

, divided by

T-s

MASE = 0.6. Did your model beat the naive in the holdout?

You cannot tell: the naive’s holdout errors never enter MASE. MASE $<1$ does not mean you beat the naive; $U_2<1$ does.

Helpers theil_u2(actual, forecast, benchmark) and mase(actual, forecast, train, m = 1): the helper’s m is the seasonal lag

s

. accuracy() differs: a seasonal MASE scale by default, and a different “Theil’s U”.

[About 4 minutes. Socratic predict-before-show: ask the question, take answers, then click for the answer.] Absolute scores need a yardstick. RMSE and MAE are in the units of y and MAPE is a percentage; an RMSE of a tenth of a point sounds good, but it depends entirely on how hard the series is to forecast. A model is worth building only if it beats something trivial, so we need relative scores. Two are in common use and they answer different questions; conflating them is the most common error in this part (Common Pitfall 1). U_2 is a race: the model’s RMSE over locf’s RMSE, on exactly the same N holdout dates. Below 1 it beat the naive there; 1 is a tie; above 1 it lost. MASE is a yardstick, not a race: the holdout MAE divided by the MAE the naive forecast made one step ahead across the training sample (Hyndman and Koehler 2006). The yardstick was measured on a different stretch of data. MASE below 1 says your holdout errors were smaller than the naive’s typical training miss. It does not say you beat the naive in the holdout: the naive’s holdout errors never enter the formula. If the holdout was an unusually calm stretch, everything, including the naive, has MASE below 1. That is the answer to the prompt. Seasonal MASE uses the seasonal naive |y_t - y_{t-s}| over T - s in the denominator. Why not just MAPE: it divides by y_t, so it is undefined at zero and explodes near it; MASE scales by a training error instead. The helpers: theil_u2() and mase(…, m = 1) in helpers/forecast_utils.R. Notation flag: the helper’s m is the seasonal lag, written s in the math, because m is the forecaster index today. accuracy() traps, one line each: for a seasonal ts it scales MASE by the seasonal naive (m = 12) by default; UNRATE is seasonally adjusted at the source, so the non-seasonal scale is the honest one. Its “Theil’s U” works with relative changes and drops the first date: a cousin of U_2, not the same number. Theil’s U_1 is a Deeper Dive; leave it. [Sources] notes 7.3 “Concept: Two Different Questions”, Math (the two formulas, reading U_2, reading MASE, seasonal MASE, why not MAPE, the helper’s m), “Other Theil Statistics, and accuracy() Defaults”; Common Pitfalls 1. Merges the scaffold’s two U_2/MASE stubs (u2-mase-concept, u2-mase-math) into one slide, per the timing brief.

The scorecard on UNRATE: small errors, no win

metrics_table <- function(actual, Fmat, train) {
  bench <- Fmat[, "locf"]
  t(apply(Fmat, 2, function(f) c(
    RMSE = rmse(actual, f),
    MAE  = mae(actual, f),
    MAPE = mape(actual, f),
    MASE = mase(actual, f, train),
    U2   = theil_u2(actual, f, bench))))
}
metrics_table(y_h, F_h, train_un)   # 2019

Every credible track has MASE below 1. Did any of them beat locf?

Track	RMSE	MAE	MAPE	MASE	$U_2$
`arima_1`	0.097	0.083	2.258	0.642	1.011
`arima_2`	0.097	0.079	2.127	0.607	1.016
`arima_3`	0.096	0.080	2.168	0.618	0.999
`locf`	0.096	0.075	2.030	0.580	1.000
`avg_k`	0.108	0.092	2.478	0.709	1.124
`meanf`	2.315	2.312	63.103	17.875	24.182

No meaningful win: arima_3’s 0.999 is a tie with locf. MASE 0.58 to 0.71, but

U_2

0.999 to 1.124 (locf is 1 by construction). 2019 was easy, for everyone, the naive included.

meanf: an I(1) series has no

\mu

to converge to (Module 2).

[About 5 minutes. LIVE CODING beat: type metrics_table() and run it on the UNRATE holdout; the table is the output, already on the slide. Ask the question before clicking. Fallback: have F_h, y_h and train_un pre-run, in case the live coding runs long.] One function builds the scorecard: a row per track, locf always the U_2 benchmark. rmse(), mae(), mape(), mase() and theil_u2() are the course helpers in helpers/forecast_utils.R; metrics_table() itself is the notes’ inline function. Read the table in the order of the prompt. Every credible track (the three ARIMAs, locf, avg_k) has MASE between 0.58 and 0.71, comfortably below 1, while U_2 for the four non-locf tracks sits between 0.999 and 1.124 (locf is 1 by construction). In words: small errors by the training yardstick, and no better than locf on these dates. The closest, arima_3 at U_2 = 0.999, is a tie for practical purposes, not a win. 2019 was a calm glide in unemployment, so every reasonable forecast had small errors, including the naive. MASE below 1 is a statement about 2019 being easy, not about the models being good. Close tracks, two losses: locf has the smallest MAE (0.075), while on RMSE it and arima_3 are essentially tied. Module 6’s lesson: the loss picks the question, and with gaps this small neither metric orders the close tracks with any confidence. That is the bridge to DM. meanf: RMSE 2.32, MASE 17.9, U_2 24.2. Module 2 callback: the order of integration decides which baselines are even credible. An expanding mean converges to mu for a stationary series; an I(1) series has no mu. It returns in 7.8. If asked about the simulation (one example on the deck; this one stays in the notes): on the course DGP the correctly specified model, arima_1, wins this realization on every metric (U_2 = 0.86), and arima_2 shows the disagreement the other way round: U_2 = 0.93 (it beat locf on those dates) with MASE = 1.00 (its miss was about the naive’s typical training miss). Both true; different questions. [Sources] notes 7.3 “Simulate: The Metrics Table on the Course DGP” (chunk metrics-fn; the three passes) and “Real Data: The Metrics Table on UNRATE” (chunk metrics-unrate; the MASE/U_2 contrast, meanf, the Module 2 callback); numbers from cache/m7_results.rds (met_un, met_sim). Merges the scaffold’s metrics-sim stub into this slide.

Scoring is not testing: the taste test

A metrics table is a ranking, not a verdict. Is the gap bigger than what the noise in twelve forecast errors would produce by chance?

Two coffee shops, one sip eachno verdict

Twelve mornings, one sip from eacha tally of the difference

Is the average difference zero?a paired

t

-test

Shop A won one sip. How many paired sips before you would bet on A?

The Diebold-Mariano test is a $t$ -test on the loss differentials.

[About 2 minutes. Concept step for section 7.4; the punchline is a fragment, after the prompt.] A metrics table tells you which number is smaller. It does not tell you whether the gap is larger than the noise in twelve forecast errors would produce by chance. Module 6’s question was exactly this: two RMSEs that differ in the second decimal place, is that a real difference? The taste test. Two coffee shops, one sip each, and you prefer shop A. Is A better? One sip is not a verdict. Twelve paired sips, one from each shop each morning, with a running tally of how much better A tasted each day: now you can ask whether the average difference is distinguishable from zero. That is a t-test on the paired differences. Take answers to the prompt; the honest answer is “it depends on how much the daily difference bounces around”, which is exactly the next slide’s denominator. Then click: that is all Diebold-Mariano is. [Sources] notes 7.4 “Concept: Scoring Is Not Testing” (the taste-test analogy).

Diebold-Mariano: a $t$ -test on loss differentials

$d_t=g(e_{1,t})-g(e_{2,t}),\qquad H_0\!:\ \mathbb{E}[d_t]=0$

$\bar d=\frac{1}{N}\sum_t d_t,\qquad \text{DM}=\frac{\bar d}{\hat\sigma_{\bar d}}\ \stackrel{a}{\sim}\ N(0,1)$

\hat\sigma_{\bar d}

: a standard error that allows for serial correlation in $d_t$ . At

h=1

, dm.test() uses lag 0 only,

\hat\sigma_d/\sqrt N

(after the small-sample correction): it assumes $d_t$ is serially uncorrelated, a condition on

d_t

, not merely on the errors.

Out-of-sample errors onlynever residuals(fit)

Same dates, same lossboth on the same

N

The sign:

\text{DM}<0

forecast 1 had the smaller loss

On the course simulation, the correctly specified model against the naive:

es_h <- ys_h - Fs_h        # N x M holdout errors
dm.test(es_h[, "arima_1"], es_h[, "locf"],
        h = 1, power = 2)  # power 2: squared
d_sim <- es_h[, "arima_1"]^2 - es_h[, "locf"]^2
t.test(d_sim)$statistic

dm.test()DM = −0.90, p = 0.390

t.test(d_sim)t = −0.90

The same statistic: dm.test() applies the Harvey-Leybourne-Newbold correction and a

t_{N-1}

reference. Named, not derived.

[About 6 minutes. Math step, then the worked example; present and explain.] Two forecasts, out-of-sample errors e_{1,t} and e_{2,t} over N dates, a loss g chosen in advance (Module 6: g(e) = e^2 or |e|; g, not L, because L is the lag operator). The loss differential d_t is the loss of forecast 1 minus the loss of forecast 2. If they are equally accurate it is zero on average: H_0: E[d_t] = 0. DM is the mean differential over its standard error, approximately standard normal under the null. The clause that matters: sigma-hat_{d-bar}, the standard error of the mean, must allow for serial correlation in d_t, because loss differentials are a time series. With h = 1, dm.test() uses the lag-0 variance only, so it becomes the familiar sigma-hat_d over root N, after the small-sample correction (sigma-hat_d the sample SD of d_t). That version does not allow for serial correlation: it assumes d_t is serially uncorrelated. Say clearly that this is a condition on the loss differential, not merely on the forecast errors. Optimal one-step forecasts have serially uncorrelated errors, which makes the assumption plausible, but it is not guaranteed: d_t = g(e_{1,t}) - g(e_{2,t}) is a nonlinear function of two error series and can be autocorrelated even when neither error is, and a misspecified competitor makes that more likely. So h = 1 assumes serially uncorrelated d_t; whether that holds for the UNRATE tests is the section 7.8 check (optional; one line of results on the unrate-dm slide). Why multi-step errors are autocorrelated, and the long-run variance, is the Technical Note; leave it there. If asked how to check (7.4 Technical Note): look at the ACF of d_t, then re-run with a variance that allows for serial correlation (h = 2 or more, varestimator = “bartlett”) and see whether the verdict moves. Three rules. Out-of-sample errors only: in-sample residuals reward overfitting. Same dates, same loss. The sign tells you who won: DM below 0 means d-bar below 0, forecast 1 had the smaller average loss. In R: dm.test(e1, e2, h = 1, power = 2); h is the horizon, always 1 today; power = 2 squared, power = 1 absolute. The worked example, the course simulation’s twelve holdout dates. Read the output line by line: DM = -0.90, negative, so arima_1 had the smaller squared errors; the horizon and power echo your choices; p = 0.390; two-sided. And the one-sample t-test on d_t returns the same statistic, -0.90: at h = 1, after the small-sample correction, DM is literally a t-test on the loss differentials. Harvey, Leybourne and Newbold (1997): name it, do not derive it (Enders eq. 2.62 and the text after it; Enders cites a different HLN paper, 1998). And the hook for the next slide: arima_1 is the correctly specified model and won on all five metrics, yet at N = 12 the test does not reject. Bad luck, or typical? [Sources] notes 7.4 Math (d_t, H_0, DM, the serial-correlation clause, “What h = 1 assumes”, three rules, dm.test syntax, HLN) and “Worked Example: DM Really Is a t-Test” (chunk dm-sim, verbatim); notation dictionary (sigma-hat_{d-bar} vs sigma-hat_d, g). Numbers: cache/m7_results.rds (dm_sim). The 7.4 Technical Note (long-run variance) and the Clark-West Deeper Dive stay off the deck; the 7.8 Technical Note’s dependence check appears, as one line of results, on dm-unrate and unrate-dm.

Shrink the sample: DM power from $N=12$ to $120$

arima_1, correctly specified and estimated, has an expected-loss advantage over locf and arima_2 on average. How often does DM detect it?

R=200

replications: a fresh series, the ARIMAs estimated once on 240 months, one-step tracks over the next 120, DM on the first

N

dates.

N=12

, what share of replications reject in favour of arima_1: about 90%, 50%, or 20%?

Rejecting for arima_1	N=12	N=24	N=60	N=120
arima_1 vs locf	12.5%	40.5%	76.0%	95.5%
arima_1 vs arima_2	10.0%	29.0%	60.5%	86.0%

Mainly power: DM rejects when the gap is large relative to the variability of

d_t

. “Failed to reject” never means “equally good”.

Caveat: these are rates of the

h=1

test, calibrated only if

d_t

is serially uncorrelated. Under a calibrated null, 2.5% would reject in this direction (dashed line).

[About 5 minutes. The assumption-break for Part A: make the sample tiny. Predict-before-show: take a show of hands on 90 / 50 / 20, then click once (table and figure appear together).] The DGP was built so that the correctly specified model, arima_1 with its coefficients estimated from 240 observations, has a real expected-loss advantage over locf and over arima_2. Estimated coefficients carry error, so that advantage is not guaranteed in every replication, but it is there on average. The question is how often DM detects it with N evaluation dates. Monte Carlo, seed 6377: R = 200 replications; each draws a fresh series, estimates the three ARIMAs once on 240 training months, builds one-step tracks over the next 120 months and runs DM on the first N in {12, 24, 60, 120}. The test is two-sided; the table counts only rejections in the right direction, arima_1 the winner. So the dashed reference line is at 2.5%, not 5%: under a calibrated null, half of the two-sided 5% rejections would go the other way. The answer: at N = 12, the PS3 holdout, DM rejects in favour of arima_1 in 12.5% of replications against locf and 10.0% against arima_2: lower than even the lowest option offered. At N = 120, 95.5% and 86.0%. The correctly specified model has the advantage on average; with twelve dates the test usually cannot detect it. The main lesson is power, with a calibration caveat. DM rejects when the loss gap is large relative to the variability of d_t; at N = 12 a real but moderate gap is usually invisible. The caveat: these are rejection rates of the h = 1 test, whose 5% level is calibrated only if d_t is serially uncorrelated. locf and arima_2 leave structure in their errors that arima_1 captures, so that is an assumption about these d_t, not a guarantee; the table measures power to the extent the assumption holds. A large gap still shows up at N = 12: section 7.8 produces exactly that, on UNRATE (seed it, do not spoil it). So “DM failed to reject” never means “the forecasts are equally good”; it means “these dates cannot separate them” (Common Pitfall 2). If asked about nesting: arima_2 is nested in arima_1, and here the large model is the correctly specified one, the alternative rather than the null. Clark-West is designed to address standard DM’s size and power problems in nested comparisons like this one; we have not run it in the simulation, so the table says nothing about how much it would change. Deeper Dive, not assessable. The loop is 200 x 4 x 2 DM tests on 200 simulated series: not live. The cached result is the notes’ own run. [Sources] notes 7.4 “Simulate - The Assumption-Break: Make the Sample Tiny” (chunks dm-power, plot-dm-power; seed 6377) and “The main lesson is power - with a calibration caveat”; the Clark-West Deeper Dive (nested comparison, not run); Common Pitfalls 2. Numbers: cache/m7_results.rds (power_tab), checked against the rendered notes by figures/make_figures.R. Figure: images/dm_power.png, from the notes’ pw data, with the reference line at the 2.5% directional rate, as in the notes’ plot-dm-power. Folds the scaffold’s dm-predict stub into this slide’s fragments.

Module 6’s question, answered

Module 6: on 2018-2019 the naive random walk beat the Module 5 SARIMA by a few hundredths. Does that margin mean anything?

Asked in the PS3 one-step design, on the 2019 holdout:

e_h <- y_h - F_h          # N x M holdout errors
dm.test(e_h[, "arima_1"], e_h[, "locf"],
        h = 1, power = 2)

arima_1 vs locfDM = 0.13, p = 0.90

Bartlett variance, h = 2:4p = 0.88 to 0.89

No: the test does not reject on these twelve months.

What it does not claim.

It does not test Module 6’s 24-step path. Errors from one origin share the early shocks and their spread grows with

h

: not the sequence DM assumes.

It does not say the two are equally good. Twelve dates rarely separate even a correctly specified estimated model from a misspecified one.

It is not assumption-free:

h=1

assumes

d_t

is serially uncorrelated. A variance that allows for serial correlation gives the same verdict (a low-power check).

[About 4 minutes. Real-data step for 7.4; closes Module 6’s open question and ends Part A.] Module 6’s result, in one line: forecasting 2018-2019 from December 2017, the Module 5 SARIMA lost to the naive random walk by a few hundredths, and Module 6 ended on a question for a statistician: does that margin mean anything? Now ask it properly. arima_1 (the Module 5 model) against locf on the 2019 holdout, PS3 one-step design. DM = 0.13, p = 0.90. Asked in this design, on 2019, the answer is no: the test does not reject on these twelve months. The verdict rests on the h = 1 assumption that these d_t are serially uncorrelated (the dm-math slide). Section 7.8’s Technical Note checks it: the lag-1 autocorrelation of this d_t is -0.33, not detectably different from zero (Ljung-Box p = 0.58, though twelve points give that check little power), and re-running with a Bartlett variance that allows for serial correlation (dm.test lag argument h = 2:4) gives p = 0.88 to 0.89. The verdict is robust to that check, not proven by it. Be careful about what it does and does not claim. It does not test Module 6’s 24-step path: that was one origin, many horizons, and DM does not apply to it. Errors made from a single origin at horizons 1 to 24 share the same early shocks, so they are serially correlated, and their spread grows with the horizon; a DM p-value on them cannot be trusted (the Technical Note shows even the MA(h - 1) correction does not cover that case). The valid design is the one used here: rolling origins, one horizon. And it does not say the two forecasts are equally good: the power simulation just showed twelve dates rarely separate even a correctly specified estimated model from a misspecified one. Optional, if asked: locf is ARIMA(0,1,0), nested in arima_1, so this “no” comes from a low-power test at N = 12, applied to a nested pair for which standard DM is not the purpose-built tool: weak evidence of equal accuracy (Clark-West Deeper Dive). And the tsCV() Deeper Dive with ten years of re-estimated origins gives a different verdict, with more dates and a different scheme. Both stay in the notes. End of Part A. Transition line: on UNRATE the credible tracks are descriptively close, differing by hundredths, and the one comparison we tested does not reject. So why pick one at all? [Sources] notes 7.4 “Real Data: Closing Module 6’s Open Question” (chunk dm-unrate; the h = 1 parenthesis; the two “does not claim” paragraphs), the Clark-West Deeper Dive’s last sentences; 7.8 “Technical Note - Checking the h = 1 Assumption on d_t” (chunk dm-dependence); 7.1 “The Module 6 Handoff” (the one-line Module 6 result); Module 6 notes 6.8. Numbers: cache/m7_results.rds (dm_un, dmdep_acf, dmdep_sens).

If you can’t pick a winner, why pick?

On UNRATE the five credible tracks are descriptively close, holdout RMSE 0.096 to 0.108, and the one tested comparison, arima_1 vs locf, does not reject. Picking the winner of a twelve-month race is mostly picking noise.

Each forecast carries some information, and the loser may carry information the winner lacks.

The portfolio analogy. Forecasts are assets; forecast errors are their risks. Imperfectly correlated errors diversify, as imperfectly correlated returns do.

Bates and Granger (1969)weight by accuracy

Granger and Ramanathan (1984)weights as a regression

Stock and Watson (2004)name the puzzle

Did the crowd mean beat most of you? Why might it?

[About 3 minutes. Part B, Combinations, starts here (why-combine to unrate-reading, about 67 minutes).] Section 7.4 left us in an awkward spot. On UNRATE the five credible tracks are descriptively close: their RMSEs differ by hundredths (holdout RMSE 0.096 to 0.108), and the one comparison we tested, arima_1 against locf, does not reject. (We did not test every pair; say “descriptively close”, not “DM cannot separate them”.) Picking the winner of a twelve-month race is mostly picking noise. The better question: if you cannot pick a winner, why pick at all? Each forecast carries some information, and the loser may carry information the winner lacks. The portfolio analogy: forecasts are assets, errors are their risks. A portfolio of imperfectly correlated assets has lower variance than its components, because some losses offset others; a combination of imperfectly correlated forecasts has smaller error variance for the same reason. History in one breath: Bates and Granger (1969) weighted forecasts by their accuracy; Granger and Ramanathan (1984) recast the weights as a regression; Stock and Watson (2004) named the regularity that drives Part B, the forecast combination puzzle. (The notes add Diebold 1989 on encompassing and Diebold-Pauly 1990 on shrinkage.) The prompt is the radar callback: check first whether it did, then why it might: your errors were imperfectly correlated, so averaging cancelled some of them. The next two slides make that formal. [Sources] notes 7.5 “Concept: From Selection to Combination” (the awkward spot, the portfolio analogy, the citation line). Numbers: cache/m7_results.rds (met_un RMSE, credible tracks).

A combined forecast is a weighted sum of forecasts

$\tilde y^{(c)}_{t\mid t-1}=w_0+\sum_{m=1}^{M}w_m\,\tilde y^{(m)}_{t\mid t-1}=w_0+\mathbf{w}'\tilde{\mathbf{y}}_t$

\tilde{\mathbf{y}}_t

the

M\times1

vector of the

M

forecasts of

y_t

w_m

\mathbf{w}

the weight on forecaster

m

; the

M\times1

weight vector

w_0

an intercept: zero for every method but one

e_{c,t}=y_t-\tilde y^{(c)}_{t\mid t-1}

the combined error;

c

is a label, not an index

Hats, stars and nothing, the DGP-versus-estimation rule applied to weights:

w

generic, or a fixed rule (

1/M

)

\hat w

estimated on the weights window

w^*

the population optimum

Weights are always

w

, never

\alpha

\beta

(taken). Enders writes

\alpha_0,\ldots,\alpha_n

[About 3 minutes. Built in full: the combo2 widget’s “w” bubble links here.] The object Part B is about. A combination takes the M forecasts of the same y_t and forms an intercept plus a weighted sum. The intercept w_0 is zero for every method except Granger-Ramanathan (Method 3). The notation rule, said once: the cardinal DGP-versus-estimation rule applies to weights. The optimal population weights, built from true error variances and covariances, are starred: w*. Weights estimated from data on the weights window carry hats: inverse-MSE, Granger-Ramanathan, constrained least squares. Unhatted w is a generic weight or a fixed rule; the simple average’s 1/M is fixed, not estimated, so it carries no hat. Estimated weights are applied to the holdout as w0-hat plus w-hat-prime y-tilde_t. Weights are w because alpha (the master-equation intercept) and beta (covariate coefficients) are taken. Enders writes the Granger-Ramanathan regression with alpha_0, …, alpha_n: flag it for anyone reading Enders section 13. The radar callback: the crowd mean was this formula with w_0 = 0 and every w_m = 1/M. [Sources] notes 7.5 (the combined-forecast definition, the hat/star convention paragraph), 7.6 opening paragraph (applied to the holdout as w0-hat + w-hat’ y-tilde_t); notation dictionary, Forecasting Notation (y-tilde^(c), w_0, w_m and w, the hat convention, e_{c,t}).

Two forecasts: the combined error variance

Weights

w

and

1-w

sum to one, so the combined error is the same mix of the two errors:

$e_{c,t}=w\,e_{1,t}+(1-w)\,e_{2,t}$

$\begin{aligned}\text{Var(}e_{c,t}\text{)}&=w^2\,\text{Var(}e_{1,t}\text{)}+(1-w)^2\,\text{Var(}e_{2,t}\text{)}\\&\quad+2w(1-w)\,\text{Cov(}e_{1,t},e_{2,t}\text{)}\end{aligned}$

Minimise over

w

(Enders 2.73):

$w^*=\frac{\text{Var(}e_{2,t}\text{)}-\text{Cov(}e_{1,t},e_{2,t}\text{)}}{\text{Var(}e_{1,t}\text{)}+\text{Var(}e_{2,t}\text{)}-2\,\text{Cov(}e_{1,t},e_{2,t}\text{)}}$

With zero covariance, the inverse-variance weight:

$w^*=\frac{1/\text{Var(}e_{1,t}\text{)}}{1/\text{Var(}e_{1,t}\text{)}+1/\text{Var(}e_{2,t}\text{)}}$

Error correlation (two indices: not the lag-

k

\rho_k

$\rho_{1,2}=\frac{\text{Cov(}e_{1,t},e_{2,t}\text{)}}{\sqrt{\text{Var(}e_{1,t}\text{)}\,\text{Var(}e_{2,t}\text{)}}}$

Variance ratio

\text{Var(}e_{2,t}\text{)}/\text{Var(}e_{1,t}\text{)}\ge1

: forecast 1 is the better one.

Diversification: equal variances

V

, zero covariance:

w^*=0.5

and

\text{Var(}e_{c,t}\text{)}=0.5V

. Averaging halves the error variance.

Correlation eats the gain. Among weights in

[0,1]

, the best mix beats forecast 1 alone until

\text{Cov(}e_{1,t},e_{2,t}\text{)}=\text{Var(}e_{1,t}\text{)}

, that is,

\rho_{1,2}=\sqrt{\text{Var(}e_{1,t}\text{)}/\text{Var(}e_{2,t}\text{)}}

Past that threshold, $w^*>1$ : the formula shorts the worse forecast. Unconstrained regression weights do exactly this (Method 3).

[About 6 minutes. Built in full: the combo2 widget’s rho_{1,2}, var_ratio, varec, wstar, inv_var, threshold and short bubbles link here. Present and explain; the derivative is one line in the notes, skip it unless asked.] Because the weights sum to one, y_t = w y_t + (1 - w) y_t, so the combined error e_{c,t} is the same combination of the individual errors (c labels the combination; it is not a forecaster index). Its variance is the portfolio formula, Enders 2.71 to 2.72. Minimise over w: Enders 2.73. With zero covariance, divide top and bottom by the product of the variances: the weight on each forecast is proportional to the inverse of its error variance (the form Enders states after 2.73, unnumbered). The less accurate forecast still gets weight, just less: the motivation for Method 2. Three readings. (1) Diversification: equal variances V, zero covariance, w* = 0.5 and the combined variance is 0.25V + 0.25V = 0.5V. (Enders says “25% of” either forecast’s variance; each forecast contributes 0.25V, the total is 0.5V.) (2) With zero covariance, inverse-variance weights. (3) Correlation eats the gain: among weights in [0, 1] the best mix’s advantage over the better forecast shrinks as the error correlation rises and vanishes once Cov(e1, e2) >= Var(e1). Divide by sqrt(Var(e1) Var(e2)) for the threshold in rho_{1,2}: sqrt(Var(e1)/Var(e2)); at the notes’ figure values (Var 1 and 1.5) that is about 0.82. Beyond it the unconstrained formula gives w* > 1: it shorts the worse forecast. Hold that thought for Granger-Ramanathan. Notation, said once: Var and Cov as Enders writes them, not sigma_1^2 and sigma_12, because sigma^2 is reserved for the innovation variance. rho_{1,2} has two indices separated by a comma, to keep it apart from the lag-k autocorrelation rho_k: written without the comma, rho_12 would read as the lag-12 autocorrelation. Convention flag: Enders writes the error as y - f in the derivation and f - y in the spread example; the sign does not change any variance or weight. The next slide is the widget: the notes’ static four-curve figure (Var 1, Var 1.5, rho in {-0.5, 0, 0.5, 0.9}) made draggable. Four-step: break shared with section 7.7, real data with section 7.8 (notes section 7.5 four-step note). [Sources] notes 7.5 “Math: Two-Forecast Portfolio Algebra” (the variance, w*, the zero-covariance form, the three readings, the rho_{1,2} definition, the “Var and Cov, not sigma” note, the convention flag); notation dictionary (rho_{1,2} and its collision warning, V as a local shorthand, e_{c,t}). The threshold in rho_{1,2} form is stated in notes 7.5, right after Cov >= Var(e1); the combo2 widget uses the same expression for its “rho_{1,2} to threshold” button.

Drag the weight, the correlation, the variance ratio

makeCombo2({refs: {
  w:          {label: "Slide: A combined forecast is a weighted sum of forecasts", href: "#/comb-weights"},
  wstar:      {label: "Slide: Two forecasts: the combined error variance",        href: "#/two-forecast-algebra"},
  rho12:      {label: "Slide: Two forecasts: the combined error variance",        href: "#/two-forecast-algebra"},
  var_ratio:  {label: "Slide: Two forecasts: the combined error variance",        href: "#/two-forecast-algebra"},
  varec:      {label: "Slide: Two forecasts: the combined error variance",        href: "#/two-forecast-algebra"},
  simple_avg: {label: "Slide: Method 1: the simple average",                      href: "#/method-simple"},
  inv_var:    {label: "Slide: Two forecasts: the combined error variance",        href: "#/two-forecast-algebra"},
  threshold:  {label: "Slide: Two forecasts: the combined error variance",        href: "#/two-forecast-algebra"},
  short:      {label: "Slide: Two forecasts: the combined error variance",        href: "#/two-forecast-algebra"}
}})

[Widget makeCombo2, widgets/combo2.qmd. About 3 minutes. Closed form, nothing simulated. Mount: Var(e1) = 1, Var(e2) = 1.5 (the notes’ section 7.5 figure), rho{1,2} = 0.5, w = 0.5.] (a) At the mount, predict: does the simple average (w = 0.5) beat forecast 1 alone? Read the diamond against the dashed line: Var(e_c) at w = 0.5 is 0.931, below 1. Then snap w to w: 0.696, Var 0.882, 11.8% below forecast 1. (b) Drag rho_{1,2} up and watch the green region (the weights that beat forecast 1) shrink. Snap rho_{1,2} to the threshold, about 0.816 at ratio 1.5: the region vanishes, no gain. (c) Push past it to 0.9 (the “section 7.5 figure, rho_{1,2} = 0.9” preset): w > 1, the regime message says it shorts the worse forecast. That is what Granger-Ramanathan’s unrestricted weights will do on twelve observations. (d) Finish on “negative rho_{1,2}”: the diversification ceiling. Errors that move in opposite directions cancel hardest. Every control carries an info bubble; its “Defined on” link jumps to the slide that defines the term (Esc closes a bubble; it is kept from reveal.js). [Sources] slides/module_07/widgets/_combo2.qmd header and demo_combo2.qmd (default-instance numbers, in-class script); notes 7.5.

Method 1: the simple average

$w_0=0,\qquad w_m=\frac{1}{M}$

Nothing is estimated, so the weights carry no hat. That is the whole point: the zero-estimation-error benchmark that every other method has to beat.

It is the radar’s crowd mean: every forecaster gets one vote,

1/M

of the result.

Every method fits on the weights window (

T_w=12

actuals and the

T_w\times M

forecasts of them) and is then applied, unchanged, to the holdout forecasts.

comb_simple <- function(y, F) {
  M <- ncol(F)   # y unused: common interface
  list(intercept = 0,
       weights = setNames(rep(1 / M, M),
                          colnames(F)))
}

[About 2 minutes. Built in full: both widgets’ “simple average” bubbles link here.] Method 1 is the benchmark. w_0 = 0 and each of the M weights is 1/M; with M = 6, one sixth each. Nothing is estimated from the window, which is why the weights carry no hat: the simple average is a fixed rule, and its estimation error is exactly zero. Every comparison in section 7.7 is relative to this method. The common contract all four methods share: input the weights-window actuals y and the T_w x M matrix of forecasts F; output list(intercept, weights). F is only an R argument name; in the math the forecast matrix gets no symbol (F is the Granger F-test in Module 9). Canonical copies in helpers/combination_utils.R, which Problem Set 3 can source. Radar callback: the crowd mean, one vote each, including the one absurd guess. Seed the caveat again: a broken member gets a full vote. [Sources] notes 7.6 opening paragraph (the weights window, the holdout, the common contract, F as an argument name only) and “Method 1 - Simple Average” (chunk comb-simple, verbatim); notation dictionary (the simple average’s 1/M is fixed and carries no hat).

Method 2: inverse-MSE weights

$\hat w_m=\frac{1/\text{MSE}_m}{\sum_{j=1}^{M}1/\text{MSE}_j},\qquad \text{MSE}_m=\frac{1}{T_w}\sum_{t\in\text{weights}}e_{m,t}^2,\qquad \hat w_0=0$

The

M

-forecast version of the zero-covariance

w^*

(Bates and Granger 1969; Enders 2.74).

$M$ numbers estimated, each from its own column. Non-negative and summing to one by construction.

MSE, not the error variance, so a biased forecaster pays for its bias too.

comb_inv_mse <- function(y, F) {
  mse <- colMeans((y - F)^2)
  list(intercept = 0,
       weights = (1 / mse) / sum(1 / mse))
}

[About 3 minutes. Built in full: the combination-puzzle widget’s “inverse-MSE” bubble links here.] Take the zero-covariance optimal weight from the two-forecast algebra and extend it to M forecasts: each forecaster’s weight is proportional to the inverse of its weights-window MSE. Bates and Granger (1969); Enders eq. 2.74. How much is estimated: M numbers, but each from its own column only, never jointly, so the estimation error is low. The weights are non-negative and sum to one by construction, with no intercept. The notes use MSE rather than the error variance on purpose: a forecaster that is consistently off gets penalised for that bias, not just its spread. Hats: these weights are estimated from the window, so w-hat. The helper version (helpers/combination_utils.R) also stops on a zero-MSE track (“a leak, not a forecast”); the slide shows the two-line core. [Sources] notes 7.6 “Method 2 - Inverse-MSE Weights” (formula, the one-column-at-a-time point, MSE vs variance; chunk comb-inv-mse, guard line omitted on the slide).

Method 3: Granger-Ramanathan

$y_t=w_0+\sum_{m=1}^{M}w_m\,\tilde y^{(m)}_{t\mid t-1}+v_t,\qquad t\in\text{weights window}$

Unconstrained least squares with an intercept. OLS on the weights window delivers

\hat w_0,\hat w_1,\ldots,\hat w_M

The intercept absorbs any bias; weights can be negative and need not sum to one.

With

M=6

and

T_w=12

: seven numbers from twelve observations.

v_t

is the regression’s disturbance: not an innovation

\epsilon_t

, not a forecast error

e_{m,t}

comb_gr <- function(y, F) {
  b <- coef(lm(y ~ F))
  w <- setNames(b[-1], colnames(F))
  w[is.na(w)] <- 0   # exactly collinear track
  list(intercept = unname(b[1]), weights = w)
}

Module 4 asked the same question of extra AR lags. Do parameters that fit the sample forecast?

[About 5 minutes. Built in full: the combination-puzzle widget’s “Granger-Ramanathan” bubble links here.] Granger and Ramanathan (1984), Enders eq. 2.75: regress the actuals on the forecasts, with an intercept and no constraints, over the weights window. The estimates are w0-hat, w1-hat, …, wM-hat. The intercept absorbs systematic bias; the weights can be negative (the shorting from the two-forecast algebra) and need not sum to one. It is one lm() call. Say it out loud: M = 6 and T_w = 12, so seven numbers from twelve observations. (On the course simulation one track is exactly collinear and R returns NA for it: arima_2 is an AR(1) with a mean, so its forecast is a constant plus a multiple of locf, and with an intercept in the regression locf adds nothing. Treat the NA as zero, the same as dropping the duplicate; leave it in and every holdout forecast is NA. On UNRATE all seven are estimated.) v_t: the regression’s disturbance, its own symbol, so it collides with nothing: not epsilon_t, not a model residual e_t, not a forecast error e_{m,t}, not the stationary deviation u_t. Enders flag: he writes these weights alpha_0, …, alpha_n. Module 4 callback for the prompt: extra parameters that fit the sample are not the same as parameters that forecast. Section 7.7 makes that concrete. [Sources] notes 7.6 “Method 3 - Granger-Ramanathan (Unconstrained Least Squares)” (the regression, v_t, the NA explanation, chunk comb-gr verbatim, “seven numbers from twelve observations”); notation dictionary (v_t, w_0 is zero except for Granger-Ramanathan).

Method 4: constrained least squares

$\hat{\mathbf{w}}=\arg\min_{\mathbf{w}}\sum_{t\in\text{weights}}\big(y_t-\mathbf{w}'\tilde{\mathbf{y}}_t\big)^2\quad\text{subject to}\quad w_m\ge0,\;\;\sum_{m=1}^{M}w_m=1,\qquad \hat w_0=0$

A quadratic program. No intercept, weights in a box: non-negative and summing to one.

The constraints do model selection: several weights land exactly on zero, and the rest read as shares.

Stable, interpretable, cannot short a forecaster.

No separate intercept for bias adjustment; reweighting can reduce bias, but cannot eliminate a common same-direction bias shared by every forecast.

comb_cls <- function(y, F) {
  M    <- ncol(F)
  Amat <- cbind(rep(1, M), diag(M))
  bvec <- c(1, rep(0, M))   # sum = 1; w >= 0
  w <- quadprog::solve.QP(crossprod(F),
         crossprod(F, y), Amat, bvec,
         meq = 1)$solution
  w <- pmax(w, 0); w <- w / sum(w)
  list(intercept = 0,
       weights = setNames(w, colnames(F)))
}

[About 5 minutes. Built in full: the combination-puzzle widget’s “constrained least squares” bubble links here.] Least squares on the weights window with no intercept, every weight non-negative, the weights summing to one. That is a quadratic program: expand the sum of squares and the objective is fully described by an M x M quadratic part (crossprod(F), F’F) and an M x 1 linear part (crossprod(F, y), F’y). quadprog::solve.QP() takes those as Dmat and dvec and the constraints as Amat and bvec, read column by column; meq = 1 makes the first column an equality. Column 1 of Amat with bvec[1] = 1 is the sum-to-one constraint; columns 2 to M+1 with zeros are the non-negativity constraints. Those are R argument names only; no math symbols, because D, d, A and b already mean other things in the course. The pmax/normalise line clips round-off like -1e-17 to zero. Reading: several weights land exactly on zero, so the inequality constraints are doing model selection; the survivors read as shares (“x% arima_1, y% avg_k”). It cannot short a forecaster. And it has no separate intercept for bias adjustment: reweighting can reduce bias (lean on the less biased forecasts), but it cannot eliminate a common same-direction bias shared by every forecast, because any weighted average of forecasts that all miss high also misses high. Housekeeping if asked: solve.QP needs crossprod(F) positive definite; two identical columns break it, and the helper version stops with a message naming them. quadprog is approved for the course (decisions log Entry 14); ForecastComb is retired. [Sources] notes 7.6 “Method 4 - Constrained Least Squares” (the program, the solve.QP mapping, Amat and bvec, the model-selection reading, the no-intercept caveat; chunk comb-cls, condensed: the constraint matrices are inlined and quadprog:: is namespaced instead of library()); design_principles section 6 (quadprog).

Four methods, one dial

Method	$w_0$	$\sum w_m=1$	$w_m\ge0$	Numbers estimated from the window	Estimation risk
Simple average (SA)	0	yes	yes	0	none
Inverse-MSE	0	yes	yes	$M$ , one column at a time	low
Constrained LS (CLS)	0	yes	yes	$M$ , jointly, boxed in	moderate
Granger-Ramanathan (GR)	free	no	no	$M+1$ , jointly, unrestricted	high

The four, fit on the course simulation’s weights window:

	simple avg	inverse-MSE	constrained LS	Granger- Ramanathan
$\hat w_0$	0.00	0.00	0.00	261.70
`arima_1`	0.17	0.27	0.66	0.78
`arima_2`	0.17	0.14	0.00	1.88
`arima_3`	0.17	0.11	0.00	−2.04
`locf`	0.17	0.09	0.00	0.00
`avg_k`	0.17	0.16	0.00	−0.83
`meanf`	0.17	0.23	0.34	−50.83

Read it as a dial. The simple average estimates nothing; inverse-MSE a little; CLS more, inside a box; Granger-Ramanathan everything.

CLS finds arima_1. Granger-Ramanathan’s intercept is 262 on a series whose mean is 5.

Which end of the dial wins on data the weights have not seen?

[About 3 minutes. The summary of section 7.6, and the bridge to the assumption-break.] Read the table down the restrictions. The simple average: no intercept, sums to one, non-negative, nothing estimated, no estimation risk. Inverse-MSE: the same restrictions, M numbers estimated but each from its own column, low risk. Constrained LS: the same box, M numbers estimated jointly, moderate risk. Granger-Ramanathan: a free intercept, no constraints, M + 1 numbers estimated jointly, high risk. The weights on the left are the notes’ weights-sim chunk: the four methods fit on the course simulation’s twelve weights-window dates (seed 6376). Inverse-MSE leans toward arima_1 (0.27) but stays near one sixth everywhere. CLS puts 0.66 on arima_1 and the rest on one other track, zero elsewhere: the constraints doing model selection. Granger-Ramanathan returns an intercept of 262 and weights as large as -50.8, on a series whose mean is 5; its locf weight is 0 because the column is exactly collinear with arima_2 plus the intercept (the NA set to zero, method-gr’s notes). The dial: SA estimates nothing, inverse-MSE a little, CLS more inside a box, GR everything. The prompt is section 7.7’s question. Take a guess from the room before moving on. [Sources] notes 7.6 “Summary” (the table, verbatim; chunk weights-sim; the dial reading). Numbers: cache/m7_results.rds (ws_sim). The full-covariance Bates-Granger Technical Note and the three-regressions Deeper Dive stay off the deck.

A combination method is an estimator

Theory says Granger-Ramanathan should win. It minimises the weights-window sum of squared errors over all intercepts and weights: no other weights can fit that window better.

But the fit is in-sample, and its weights are estimates from twelve observations. Estimates carry error.

The assumption we break: that you know the weights.

The experiment is fair. The course DGP has a dominant track, arima_1, so the optimal weights are far from

1/M

A method that finds arima_1 should beat the simple average. Can it, with

T_w=12

[About 2 minutes. Concept step of the assumption-break (section 7.7).] Theory says Granger-Ramanathan should win. It minimises the weights-window sum of squared errors over all intercepts and weights, so no other set of weights, equal, inverse-MSE or constrained, can fit that window better. But “the fit” is in-sample. Its weights are estimates, made from twelve observations, and estimates carry error. The assumption we are about to break is the one buried in every optimal-weight formula: that you know the weights. Shrink the weights window and watch. Callback to the course-simulation slide: the DGP was built so the experiment is fair. If equal weights were already optimal, every estimated weight would be pure noise and the simple average would win trivially. Here one track is the correctly specified model and the optimal weights are far from 1/M, so a weighting scheme that finds arima_1 should beat the simple average. The question is whether it can find it with twelve observations. Three stages follow: one realization (Stage 1), 200 of them (Stage 2), and a sweep over the window length (Stage 3, the widget). [Sources] notes 7.7 “Concept: A Combination Method Is an Estimator” (the theory, the assumption, “The DGP was built so the experiment is fair”).

Stage 1: weights that win the window they were fit on

RMSE	Weights window	Holdout
`arima_1`	0.706	1.120
`arima_2`	0.979	1.208
`arima_3`	1.115	1.372
`locf`	1.208	1.301
`avg_k`	0.922	1.793
`meanf`	0.769	1.701
simple avg	0.877	1.284
inverse-MSE	0.804	1.281
constrained LS	0.681	1.237
Granger-Ramanathan	0.282	3.009

One realization (seed 6376),

T_w=12

: the §7.6 weights, scored on the window they were fit on and on the holdout.

Granger-Ramanathan, weights window0.282

simple average; arima_10.877; 0.706

Granger-Ramanathan, holdout3.009

Best in the window, worst on the holdout: almost as many weights as observations, fitting the noise.

The other three sit close together; CLS, which found arima_1, is the best of them.

[About 4 minutes. LIVE CODING candidate, cheap: apply_w() and the stage1 rbind on the seed-6376 realization.] Score the section 7.6 weights on both windows. apply_w(w, F) is w0 plus F times the weights; the table is the notes’ stage1 chunk: RMSE of each single track and each combination, on the weights window and on the holdout. Granger-Ramanathan has by far the best weights-window RMSE, 0.282, against 0.877 for the simple average and 0.706 for arima_1, the correctly specified model. On the holdout it posts 3.009: the worst number in the table. The weights on the last slide explain why: an intercept of 262 and weights like -50.8 and 1.9 are not a discovery about the forecasters; they are a regression with almost as many parameters as data points, fitting the noise in twelve months. The other three methods sit close to one another (1.284, 1.281, 1.237 for SA, inverse-MSE, CLS), and CLS, which found arima_1, is the best of them; still worse than arima_1 alone (1.120). One realization could be a fluke. Stage 2 says it is not. The next slide defines the score Stage 2 uses. [Sources] notes 7.7 “Stage 1 - One Realization, T_w = 12” (chunk stage1; the reading). Numbers: cache/m7_results.rds (stage1, ws_sim).

Scoring a method: holdout RMSE relative to the simple average

$\text{relative RMSE}_{\text{method}}=\frac{\text{RMSE of the method's combination on the holdout}}{\text{RMSE of the simple average on the same holdout}}$

Below 1beat the simple average

Equal to 1a tie

Above 1lost to it

Weights fit on the weights window; scored on the 12 holdout dates the weights never saw.

Not

U_2

U_2

divides by locf’s RMSE; this ratio divides by the simple average’s. Same dates, a different benchmark.

Over $R=200$ simulated series, one ratio per replication. Three summaries:

Medianthe typical replication

90th percentilethe bad case: the right tail

Share beating the simple averagefraction of replications below 1

[About 3 minutes. Built in full: the combination-puzzle widget’s rel_rmse, median, p90 and share bubbles link here.] The scoreboard for the assumption-break. For each replication: fit the four methods on the weights window, apply the weights unchanged to the 12-date holdout, compute each combination’s holdout RMSE, and divide by the simple average’s holdout RMSE on the same dates. Below 1 means the method beat the simple average there. Across the R = 200 replications each method has 200 ratios. The notes summarise them three ways: the median (the typical replication), the 90th percentile (how bad the bad case is: the long right tail is the puzzle in one picture), and the share of replications in which the ratio is below 1. The notes’ summary function also reports the mean; the widget and the slide use the three above. The contrast to make out loud, because students will conflate them: Theil’s U_2 (section 7.3) is a model’s RMSE over locf’s RMSE; this ratio is a combination’s RMSE over the simple average’s. Both are “same dates, relative to a benchmark”; the benchmark differs. Next: predict where Granger-Ramanathan lands, then the Stage 2 box plot. [Sources] notes 7.7 “Stage 2 - Monte Carlo, R = 200” (the paragraph defining the record: “each method’s holdout RMSE relative to the simple average (below 1: beat the simple average)”; the summ() function: median, mean, p90, share_beating_SA), 7.3 (U_2 definition).

Stage 2: 200 replications at a 12-date weights window

Each replication (seed 6378): a fresh series, the ARIMAs estimated once, six tracks, four methods fit on

T_w=12

, scored on the 12-date holdout relative to the simple average.

Granger-Ramanathan beats the simple average in about 80%, 50%, or 15% of replications?

Every ratio above 5× is drawn (orange) in the shaded band: 12 of Granger-Ramanathan’s 200, the largest 11.7×; in the widget they fill the hatched lane. The summaries use every replication.

	Median	90th pct	Share < 1
inverse-MSE	0.99	1.01	74.5%
constrained LS	0.95	1.12	64.0%
Granger-Ramanathan	1.68	4.07	16.0%

The long right tail is the puzzle in one picture. Inverse-MSE is the simple average with a tilt (median 0.99); CLS beats the simple average at the median (0.95).

[About 4 minutes. Built in full: the combination-puzzle widget’s “clipped” bubble links here. Predict-before-show: take the vote on 80 / 50 / 15, then click once.] Set-up. For each of R = 200 replications: simulate a fresh series from the course DGP, estimate the three ARIMAs once, build the six tracks, fit the four combinations on the 12-date weights window, and score the 12-date holdout. The record is each method’s holdout RMSE relative to the simple average (the previous slide): below 1, it beat the simple average. Seed 6378; the loop is too slow to run live, so the result is the notes’ own run, cached. The answer: Granger-Ramanathan beats the simple average in 16.0% of replications; it lost 84.0% of the time, so Stage 1 was typical. Its median holdout RMSE is 1.68 times the simple average’s and its 90th percentile 4.07 times. The long right tail is the puzzle in one picture. The other two tell the rest of the story. Inverse-MSE is the simple average with a slight tilt (median 0.99, beats it in 74.5% of replications): it leans toward arima_1, estimating so little that it almost never does much worse. Constrained LS sits between the two on risk, a thin right tail, and in this DGP it already beats the simple average at the median (0.95), because its constraints shrink the estimation error while still letting it find the dominant track. That is evidence for the estimation-error explanation of the puzzle, not against it. The clip, because the widget uses it: the box plot is on a log scale with the full range, and the shaded band marks ratios above 5x; every ratio above 5x is drawn in orange, whether or not the box-plot rule would call it an outlier (navy points at or below 5x are the box plot’s outliers). 12 of Granger-Ramanathan’s 200 replications land there (the largest 11.7x). In the widget those ratios pile into a hatched overflow lane with their count; they are clipped only in the picture: the median, 90th percentile and share use every replication. Next: the widget sweeps T_w. Same 200 series, longer windows. [Sources] notes 7.7 “Stage 2 - Monte Carlo, R = 200” (chunks stage23, plot-stage2; seed 6378; the reading of all three methods); widgets/_cpuzzle.qmd (the clip bubble: x axis 0.5x to 5x, overflow lane). Numbers: cache/m7_results.rds (stage2, gr12_over5, gr12_max, gr12_lose), checked against the rendered notes and against the widget payload’s medians by figures/make_figures.R and widgets/make_cpuzzle_data.R. Figure: images/stage2_box.png, the notes’ s2 data, horizontal. Folds the scaffold’s stage2-predict stub into this slide’s fragment.

Stage 3: how long a window do estimated weights need?

makeCombPuzzle({refs: {
  Tw:         {label: "Slide: The PS3 design: many origins, one horizon",                   href: "#/ps3-design"},
  simple_avg: {label: "Slide: Method 1: the simple average",                                href: "#/method-simple"},
  imse:       {label: "Slide: Method 2: inverse-MSE weights",                               href: "#/method-imse"},
  gr:         {label: "Slide: Method 3: Granger-Ramanathan",                                href: "#/method-gr"},
  cls:        {label: "Slide: Method 4: constrained least squares",                         href: "#/method-cls"},
  rel_rmse:   {label: "Slide: Scoring a method: holdout RMSE relative to the simple average", href: "#/rel-rmse"},
  median:     {label: "Slide: Scoring a method: holdout RMSE relative to the simple average", href: "#/rel-rmse"},
  p90:        {label: "Slide: Scoring a method: holdout RMSE relative to the simple average", href: "#/rel-rmse"},
  share:      {label: "Slide: Scoring a method: holdout RMSE relative to the simple average", href: "#/rel-rmse"},
  clip:       {label: "Slide: Stage 2: 200 replications at a 12-date weights window",       href: "#/stage2"},
  true_model: {label: "Slide: The course simulation: one track is genuinely better",        href: "#/course-dgp"}
}})

[Widget makeCombPuzzle, widgets/_cpuzzle.qmd; the payload is inline (the notes’ own Stage 2-3 Monte Carlo, seed 6378, R = 200, via make_cpuzzle_data.R). About 3 minutes. Mount at T_w = 12.] (a) At T_w = 12, ask: which row has its cloud right of the dashed line, and how far? Granger-Ramanathan: median 1.68x the simple average, 90th percentile 4.07x, beats it in 16.0% of replications; 12 replications beyond 5x sit in the hatched overflow lane. Constrained least squares 0.95, inverse-MSE 0.99: already left of the line. (b) Predict the T_w at which Granger-Ramanathan’s median crosses 1. Step to 24 (1.12), 60 (0.98: the crossover), 120 (0.93; constrained least squares 0.90; arima_1, the correctly specified model, alone 0.90, the faint dotted line). Each dot is the same simulated series at every T_w: watch them slide. (c) The lesson: estimated weights earn their keep only with data. PS3 gives you twelve dates. Design note if asked: in the sweep the holdout moves with T_w (a deliberate departure from the PS3 split, which saves refitting the ARIMAs four times per replication). The DGP is stationary, so moving the holdout does not change the distribution of y_t. It is not completely innocuous, though: a later holdout also lengthens the history behind meanf, the expanding mean, so its forecasts are a little more precise at T_w = 120 than at T_w = 12, even on this stationary DGP. The ARIMA tracks, whose coefficients stay fixed, are unaffected. Every label carries an info bubble with a “Defined on” link back to the defining slide. [Sources] slides/module_07/widgets/_cpuzzle.qmd header and demo_cpuzzle.qmd (default-instance numbers, in-class script); notes 7.7 “Stage 3 - The Sweep” and its design note.

The forecast combination puzzle

Theory says optimally weighted combinations should dominate. In practice, the simple average is very hard to beat. Stock and Watson (2004)

Three lessons from the simulation.

1. A combination method is an estimator. Holdout performance = the population gain from its weights minus the estimation error in them.

2. The puzzle is a comparison: attainable gain against estimation cost. It bites when the gain over equal weights is small relative to the cost of estimating the weights: near-equal optimal weights, or redundant forecasts. Twelve weights-window dates make the cost high, as in most applied work, and PS3.

3. Restrictions are a dial. Moving down from

1/M

pays only when the window is long or the gain is large.

The index fund. Active managers (estimated weights) should beat the index (equal weights): they pick the winners. In practice fees eat the edge. Here the fee is estimation error.

But an index fund that holds a bankrupt stock at full weight is not a good index. Equal weights are robust to estimation error, not to a broken member of the pool.

[About 4 minutes. Merges the scaffold’s three-lessons and puzzle stubs. The “But” is a fragment: it sets up 7.8.] Stock and Watson (2004) gave it a name: theory says optimally weighted combinations should dominate; in practice the simple average is very hard to beat. The simulation shows the explanation (Timmermann 2006; Smith and Wallis 2009): the optimal weights have to be estimated, and with short windows the estimation error swamps the gain. The governing comparison is the attainable improvement over equal weights versus the cost of estimating the weights. The most important case in which the gain is small is optimal weights that are near equal, so that there is little to swamp; but that is a special case, not a necessary condition: substantially different optimal weights can still yield little gain when the forecasts are redundant, carrying much the same information (Elliott 2011, reviewed in Wang et al. 2023, section 2.6; on redundant pools, section 2.5). Enders’ spread example illustrates the near-equal case, seven forecasts with nearly identical error variances; one forecast date, so an illustration, not evidence. Three lessons. (1) A combination method is an estimator: its holdout performance is the population gain from its weights minus the estimation error in them, and the weights window is where that error is paid. (2) The puzzle is a comparison, attainable gain against estimation cost (Smith and Wallis 2009; Claeskens, Magnus, Vasnev and Wang 2016; Wang et al. 2023, section 2.6): it bites when the improvement optimal weights could deliver over equal weights is small relative to the cost of estimating them. Near-equal optimal weights are the important special case; redundant forecasts are another (Elliott 2011, reviewed in Wang et al. 2023, section 2.6; on redundant pools, section 2.5). With twelve weights-window observations the estimation cost is high, which describes most applied work, and Problem Set 3. (3) Restrictions are a dial: 1/M estimates nothing, inverse-MSE a little, CLS more inside a box, Granger-Ramanathan everything; moving down the list pays only when the window is long or the gain is large. The widget showed GR needs a window of 24 to 60 dates to cross the simple average; CLS and inverse-MSE did not. The index fund. Active managers, estimated weights, should beat the index, equal weights, in theory: they pick the winners. In practice fees eat the edge, and here the fee is estimation error. Then click: but an index fund that holds a bankrupt stock at full weight is not a good index. Radar callback: the one absurd guess got a full vote in the crowd mean. Equal weights are robust to estimation error, not to a broken member of the pool. Section 7.8 shows what that means on real data. [Sources] notes 7.7 “Three Lessons” and “The Forecast Combination Puzzle” (Stock-Watson’s statement, the explanation, the governing comparison and its near-equal and redundant cases, citing Elliott 2011, reviewed in Wang et al. 2023 section 2.6, and on redundant pools section 2.5; Enders’ spread example, the index-fund analogy and its caveat). The Looking Ahead on rolling-window weights stays off the deck.

UNRATE: weights fit on 2018

ws_un <- list(SA   = comb_simple(y_w, F_w),
              IMSE = comb_inv_mse(y_w, F_w),
              CLS  = comb_cls(y_w, F_w),
              GR   = comb_gr(y_w, F_w))
sapply(ws_un, function(w)
  c(w0 = w$intercept, w$weights))
kappa(cbind(1, F_w), exact = TRUE)

Condition number, GR design matrix41,620

Before the holdout: which of the four will the simple average beat on 2019?

	simple avg	inverse-MSE	constrained LS	Granger- Ramanathan
$\hat w_0$	0.000	0.000	0.000	−141.450
`arima_1`	0.167	0.235	0.690	1.224
`arima_2`	0.167	0.158	0.000	−5.825
`arima_3`	0.167	0.229	0.000	1.088
`locf`	0.167	0.179	0.000	5.098
`avg_k`	0.167	0.198	0.310	−2.605
`meanf`	0.167	0.001	0.000	24.803

Over twelve months the expanding mean of 696+ observations barely moves: the meanf column is almost a constant, almost collinear with the intercept.

[About 4 minutes. LIVE CODING beat: the PS3 workflow on a series that isn’t yours. The four comb_*() functions from section 7.6, applied to UNRATE’s 2018 weights window (F_w, y_w built on the six-tracks slide).] This is the complete Problem Set 3 workflow on UNRATE: N_total = 720, T = 696 (train ends December 2017), weights window 2018, holdout 2019, the six tracks from Part A. The metrics table and the Module 6 DM test are done; what remains is the combination step. Read across the rows of the weight table. Inverse-MSE gives meanf a weight of 0.001, effectively zero, and spreads the rest between 0.16 and 0.24 across the five credible tracks. CLS puts 0.69 on arima_1 and 0.31 on avg_k, zero everywhere else. Granger-Ramanathan returns an intercept of -141 and a weight of 24.8 on meanf, with large offsetting weights elsewhere. Why: over twelve months the expanding mean of more than 696 observations barely moves, so the meanf column is almost a constant, almost collinear with the intercept. The condition number of the regression’s design matrix (intercept plus six tracks; kappa(…, exact = TRUE)) is about 41,620; anything in the thousands means the least-squares solution is extremely sensitive to the data. The regression is trading an enormous intercept against an enormous meanf weight to fit twelve points. Take predictions on the prompt before moving on (the textbook answer is “all of them”; hold it). [Sources] notes 7.8 “Step 1: Fit the Four Methods on 2018” (chunk unrate-weights; the three readings). Numbers: cache/m7_results.rds (w_tab, kappa_w), checked against the rendered notes.

UNRATE: scored on 2019

2019	Holdout RMSE
`arima_1`	0.097
`arima_2`	0.097
`arima_3`	0.096
`locf`	0.096
`avg_k`	0.108
`meanf`	2.315
simple avg	0.416
inverse-MSE	0.093
constrained LS	0.094
Granger-Ramanathan	0.463

Simple average0.416

credible single tracks0.096 to 0.108

Inverse-MSE; CLS0.093; 0.094

Granger-Ramanathan0.463

Of the four combinations, the simple average beat only Granger-Ramanathan.

Inverse-MSE and CLS match or edge out the best single tracks; the simple average is roughly four times their RMSE.

[About 3 minutes. The answer to the last slide’s prompt.] apply_w() puts each method’s 2018 weights on the 2019 forecasts; the table is the holdout RMSE of the six tracks and the four combinations (notes chunk unrate-holdout). The simple average has a holdout RMSE of 0.416, roughly four times the credible single tracks’ 0.096 to 0.108. Inverse-MSE (0.093) and CLS (0.094) match or edge out the best single tracks. Granger-Ramanathan (0.463) is worse even than the simple average. So the textbook answer to the prompt, “the simple average beats all of them”, is wrong here: of the four combinations it beat only Granger-Ramanathan (among single tracks it beat only meanf). Is that a real difference, or twelve dates of noise? Test it, declared in advance. [Sources] notes 7.8 “Step 2: Score Everything on 2019” (chunk unrate-holdout; the reading). Numbers: cache/m7_results.rds (rmse_un).

UNRATE: each estimated-weight combination against the simple average

Three tests, declared in advance, all reported.

comb_h <- sapply(ws_un, function(w)   # N x 4
  w$intercept + F_h %*% w$weights)
ec_h <- y_h - comb_h      # combined errors
for (m in c("IMSE", "CLS", "GR"))
  print(dm.test(ec_h[, m], ec_h[, "SA"],
                h = 1, power = 2))

Never pick the holdout winner and then test it on the same twelve dates: that chooses the hypothesis after seeing the data.

Against the simple average	DM	$p$
Inverse-MSE	−9.1	$1.9\times10^{-6}$
CLS	−9.1	$1.8\times10^{-6}$
Granger-Ramanathan	0.87	0.40

Inverse-MSE and CLS: a clear rejection at

N=12

, the power lesson from the other side. Granger-Ramanathan: no detectable difference.

The

h=1

check, all four UNRATE tests: Ljung-Box

p

= 0.58 to 0.97 on

d_t

(little power at

N=12

); Bartlett variance, h = 2:4: no verdict changes.

[About 3 minutes. Pays off section 7.4’s forward reference: “a large gap still shows up at N = 12”.] Three tests, declared in advance and all reported: each estimated-weight combination against the simple average, squared loss, h = 1. We do not pick the holdout winner and then test it on the same twelve dates; that is choosing the hypothesis after seeing the data. Inverse-MSE against the simple average: DM = -9.1, p = 1.9e-06. CLS: DM = -9.1, p = 1.8e-06. A clear rejection at N = 12: the section 7.4 power lesson from the other side. DM rejects when the gap is large relative to the noise in d_t, and a gap this large shows up even in twelve months (Common Pitfall 3: DM does not always fail with twelve observations). Granger-Ramanathan against the simple average: p = 0.40, no detectable difference. Negative DM: forecast 1, the estimated-weight combination, had the smaller squared errors. The bottom line is the 7.8 Technical Note, in one sentence. Like every h = 1 result today, these p-values assume each d_t is serially uncorrelated. For the four UNRATE tests (these three and arima_1 against locf) the lag-1 autocorrelations of d_t run from -0.33 to 0.10, not detectably different from zero: Ljung-Box p from 0.58 to 0.97. But twelve points give that check little power, so it is consistent with uncorrelated d_t without establishing it. Re-running each test with a Bartlett (Newey-West) variance that allows for serial correlation (dm.test lag argument h = 2:4): no verdict at 5% changes (inverse-MSE and CLS stay below 2e-06; Granger-Ramanathan 0.46 to 0.56). Robust to that check, not definitive. If asked why Bartlett and not the default: the truncated “acf” sum can turn negative on twelve points, and dm.test() then warns and falls back to h = 1 (the output then reports horizon 1). The loop on the slide is a compact version of the notes’ sapply over the three methods; same tests. comb_h is the notes’ unrate-holdout chunk (apply_w() over the four weight sets), written out on the slide. [Sources] notes 7.8 “Step 3: Test Each Estimated-Weight Combination Against the Simple Average” (chunk unrate-dm-comb; the reading) and “Technical Note - Checking the h = 1 Assumption on d_t” (chunk dm-dependence; its three readings); Common Pitfalls 3. Numbers: cache/m7_results.rds (dm_comb, dmdep_acf, dmdep_sens), checked against the rendered notes by figures/make_figures.R.

The honest reading: a broken member, not a lucky estimate

Here the attainable gain is large. meanf has no business on an I(1) series, and

1/M

gives it a sixth of the vote: enough to drag the average about 0.39 points off the actuals every month.

Inverse-MSE and CLS won by down-weighting it, which is exactly what estimated weights are for: that gain pays for estimating the restricted weights.

Granger-Ramanathan lost for the Stage 1 reason, amplified by near-collinearity: the intercept and the nearly constant meanf column are almost the same regressor.

Caveats: one window, a calm glide with no turning point; the

h=1

verdicts survive a low-power dependence check, so robust, not definitive. One realization is one draw; the Monte Carlo is the distribution.

For Problem Set 3: look at your six tracks before you combine them. If one is broken on your series, say so, explain why from the series’ properties, and expect it to drag the simple average.

Equal weights are robust to estimation error. What are they not robust to?

A broken member of the pool.

[About 4 minutes. The real-data step of Part B, read honestly; ends Part B.] This is not the textbook story, and the notes would be dishonest to pretend otherwise. The simple average did not win, and the reason is not that the estimated weights got lucky. Go back to the puzzle’s governing comparison, attainable gain against estimation cost: here the attainable gain over equal weights is large enough to pay for estimating the weights. The pool contains meanf, a forecaster with no business on an I(1) series, and 1/M gives it a sixth of the vote: enough to drag the average about 0.39 percentage points off the actuals every month. Equal weights are robust to estimation error, not to a broken member. Inverse-MSE and CLS won by down-weighting meanf, which is exactly what estimated weights are for: the gain from sidelining one broken member dwarfs the estimation cost of these restricted weights. The tests rest on an assumption about d_t, and the check qualifies rather than confirms it. Every DM test here used h = 1, which assumes serially uncorrelated loss differentials; the four d_t series show no detectable autocorrelation (Ljung-Box p at least 0.58), but twelve dates give that check little power, and with a serial-correlation-robust variance no verdict changes. So the UNRATE verdicts are robust to that check, not definitive. Granger-Ramanathan lost for the Stage 1 reason, amplified by near-collinearity: seven numbers from twelve observations, two of which, the intercept and the nearly constant meanf column, are almost the same regressor. Caveats. One window. 2018-2019 was a calm glide with no turning point, so every credible track had small errors; a different window could rank the methods differently. One realization is one draw; the Monte Carlo of section 7.7 is the distribution. For Problem Set 3: look at your six tracks before you combine them. If one is broken on your series (meanf on an I(1) series is the usual suspect), say so, explain why from the series’ properties, and expect it to drag the simple average. The scorecard still reports every track and every combination; the reading is where the judgment goes. The prompt’s answer, then click: a broken member of the pool. Radar callback: the one absurd guess. If asked “could you just drop it, or use a median?”: that is the Deeper Dive on robust combinations (median, trimmed mean, the simple average of the five credible tracks); it is optional and stays off the deck. [Sources] notes 7.8 “The Honest Reading” (its bullets: the attainable gain, the d_t assumption, Granger-Ramanathan, caveats) and “What this means for Problem Set 3”; 7.7 “The Forecast Combination Puzzle” (the governing comparison); Common Pitfalls 7. Numbers: cache/m7_results.rds (meanf_drag, dmdep_acf).

Key takeaways + next time

1. Many origins, one horizon. Estimate once; roll the information set.

2. $U_2$ : did you beat the naive here? MASE: how big was the miss against the naive’s training miss?

3. A metrics table is a ranking, not a verdict. DM is a

t

-test on

d_t

(

h=1

: assumes serially uncorrelated

d_t

); twelve dates buy little power. UNRATE: five tracks descriptively close; arima_1 vs locf does not reject.

4. Combine because errors diversify. Zero covariance: weights

\propto

inverse error variance.

5. A combination method is an estimator. With

T_w=12

, estimate as little as you can.

6. Equal weights are the benchmark to beat, when the pool is sound. The puzzle: estimation cost outweighs the attainable gain.

7. Every forecast here used only the series’ own past.

Next time, Module 8: ADL models.

What if an exogenous regressor could improve the forecast?

$y_t=\alpha+\beta x_t+\sum_{j=1}^{p}\phi_j y_{t-j}+\epsilon_t$

A new slider on the mixing board:

\beta x_t

and its lags. Combinations only reshuffle the information already in the pool.

[About 3 minutes. The required closing slide.] The seven Key Takeaways from the notes, condensed. (1) Many origins, one horizon: the PS3 design estimates each model once and rolls the information set; every forecast is one step ahead; three k-period average forecasts sit beside three ARIMA tracks. (2) U_2 asks “did you beat the naive here?”; MASE asks “how big was your miss relative to the naive’s typical miss in training?” Report both; never read MASE below 1 as a win. (3) A metrics table is a ranking, not a verdict: DM is a t-test on the loss differentials d_t, with a standard error that allows for serial correlation; at h = 1, dm.test() uses the lag-0 version, which assumes d_t is serially uncorrelated. Whether it rejects depends on the gap relative to the noise in d_t, and twelve dates buy little power. On UNRATE 2019 the five credible tracks are descriptively close (RMSE), and the one track-vs-track comparison tested (arima_1 vs locf) does not reject. (4) Combine because errors diversify; with zero error covariance the optimal weight is proportional to the inverse error variance. (5) A combination method is an estimator: with twelve weights-window observations, estimate as little as you can; Granger-Ramanathan’s seven numbers win the window and lose the holdout. (6) Equal weights are the benchmark to beat when the pool is sound: the forecast combination puzzle is the cost of estimating the weights outweighing the attainable gain over equal weights, most often because the optimal weights are near equal or the forecasts are redundant. A broken pool member like meanf on an I(1) series makes the attainable gain large, so estimated weights (or robust ones, a Deeper Dive) can pay for themselves. (7) Every forecast here used only the series’ own past. Next time, the handoff in the course overview’s words: what if an exogenous regressor could improve the forecast? On UNRATE in 2019 the five credible tracks were descriptively close in RMSE, the one track-vs-track comparison tested (arima_1 vs locf, carrying the last value forward) did not reject, and the best combinations only matched the best single tracks’ RMSE. Combinations can only reshuffle the information already in the pool. Could another variable, initial unemployment-insurance claims, say, carry information no univariate track has? Module 8, ADL Models and Dynamic Multipliers, turns on a new slider, beta x_t and its lags, and asks how a change in x today works its way into y over time. The equation is a simplified form, the notation dictionary’s ARMAX row (AR terms plus beta x_t, without the trend and MA terms of the full master equation); Module 8 adds the lags of x and fixes their notation. [Sources] notes “Key Takeaways” 1-7 and “Looking Ahead - Module 8”; course_overview_15kft.md L7 forward handoff (“What if an exogenous regressor could improve the forecast?” -> L8 ADL) and the handoff map (“L7 combination puzzle -> motivates L8 add exogenous information directly via ADL”).

Interactive lab machinery

This uncounted support slide keeps both widget factories in the document. Use the RevealJS menu to return to a widget if a browser reload interrupts an interaction.

// ============================== combo2Core ===============================
// Closed-form two-forecast algebra (notes section 7.5). v1 = Var(e_1), v2 = Var(e_2).
combo2Core = {
  const cov = (v1, v2, rho) => rho * Math.sqrt(v1 * v2);
  const varC = (w, v1, v2, rho) => {
    const c = cov(v1, v2, rho);
    return w * w * v1 + (1 - w) * (1 - w) * v2 + 2 * w * (1 - w) * c;
  };
  const denom = (v1, v2, rho) => v1 + v2 - 2 * cov(v1, v2, rho);
  // Enders (2.73)
  const wStar = (v1, v2, rho) => (v2 - cov(v1, v2, rho)) / denom(v1, v2, rho);
  const varStar = (v1, v2, rho) => { const c = cov(v1, v2, rho); return (v1 * v2 - c * c) / denom(v1, v2, rho); };
  // zero-covariance (inverse-variance) weight on forecast 1
  const wInv = (v1, v2) => (1 / v1) / (1 / v1 + 1 / v2);
  // rho_{1,2} at which Cov = Var(e_1) (forecast 1 the better one): no gain inside [0, 1]
  const rhoThreshold = (v1, v2) => Math.sqrt(v1 / v2);
  // Var(e_c)(w) = Var(e_1) has roots w = 1 and w = (v2 - v1)/D; the combination beats
  // forecast 1 strictly between them (empty when they coincide).
  const beatInterval = (v1, v2, rho) => {
    const r2 = (v2 - v1) / denom(v1, v2, rho);
    return r2 <= 1 ? [r2, 1] : [1, r2];
  };
  return { cov, varC, denom, wStar, varStar, wInv, rhoThreshold, beatInterval };
}

// ============================== makeCombo2 ===============================
makeCombo2 = function (opts) {
  opts = opts || {};
  const uid = combo2Uid();
  const K = combo2Core;

  // ---- grids, state ------------------------------------------------------
  // w's top end, 2.25, clears the largest w* the rho and ratio sliders can produce: 2.101 at
  // ratio 1.9, rho_{1,2} = 0.95 (sup over the continuous ranges 2.1013 at ratio 1.908). The
  // threshold button never takes rho_{1,2} above 0.95 (thrAvail below), so that bound holds.
  const G = { w: { min: -0.5, max: 2.25, step: 0.01 },
              rho: { min: -0.95, max: 0.95, step: 0.01 },
              ratio: { min: 1, max: 6, step: 0.05 } };
  const snap = (x, g) => {
    const v = Math.min(g.max, Math.max(g.min, Number(x)));
    return Math.round((v - g.min) / g.step) * g.step + g.min;
  };
  const clean = x => Math.round(x * 1e6) / 1e6;
  const DEF = { w: clean(snap(opts.w0 == null ? 0.5 : opts.w0, G.w)),
                rho: clean(snap(opts.rho0 == null ? 0.5 : opts.rho0, G.rho)),
                ratio: clean(snap(opts.ratio0 == null ? 1.5 : opts.ratio0, G.ratio)) };
  const NOTES = opts.layout === "notes";
  const SHOW_PRESETS = opts.presets !== false;
  const V1 = 1;                                 // Var(e_{1,t}) fixed
  let w = DEF.w, rho = DEF.rho, ratio = DEF.ratio;
  const AT_TOL = 0.005;                         // half a rho_{1,2} step: "at threshold" band
  // threshold button: only where the threshold sqrt(Var(e_1)/Var(e_2)) is on the rho_{1,2}
  // slider (ratio > 1.10). At ratios 1.05 and 1.10 it is 0.976 and 0.953; allowing those let
  // a later ratio drag push w* to 2.79, off the w axis.
  const thrAvail = () => ratio > 1 && K.rhoThreshold(V1, ratio * V1) <= G.rho.max + 1e-12;
  // w* outside the w axis: unreachable from the controls (see G.w), kept as a flagged fallback
  const wsOff = x => x > G.w.max + 1e-12 || x < G.w.min - 1e-12;
  const clampW = x => Math.min(G.w.max, Math.max(G.w.min, x));

  // ---- "Defined on:" pointers (opts.refs merges over these, key by key) --------
  const REF_DEFAULTS = {
    w:          { label: "Notes §7.5 — combination weights (Concept: From Selection to Combination)" },
    wstar:      { label: "Notes §7.5 — combination weights; w* is Enders (2.73), Math: Two-Forecast Portfolio Algebra" },
    rho12:      { label: "Notes §7.5 — two-forecast algebra (the error correlation is defined above the figure)" },
    var_ratio:  { label: "Notes §7.5 — two-forecast algebra" },
    varec:      { label: "Notes §7.5 — two-forecast algebra, Enders (2.71)–(2.72)" },
    simple_avg: { label: "Notes §7.6 — Method 1, simple average" },
    inv_var:    { label: "Notes §7.5 (zero-covariance w*) and §7.6 — Method 2, inverse-MSE weights" },
    threshold:  { label: "Notes §7.5 — “Correlation eats the gain”" },
    short:      { label: "Notes §7.5 — shorting the worse forecast; §7.6 — Method 3, Granger-Ramanathan" }
  };
  const REFS = {};
  Object.keys(REF_DEFAULTS).forEach(k => { REFS[k] = Object.assign({}, REF_DEFAULTS[k], (opts.refs || {})[k] || {}); });
  Object.keys(opts.refs || {}).forEach(k => { if (!REFS[k]) REFS[k] = Object.assign({}, opts.refs[k]); });

  // ---- self-checks (notes Practice Problem 4) -------------------------------
  const problems = [];
  const near = (a, b, name) => { if (!(Math.abs(a - b) < 1e-9)) problems.push(`${name}: got ${a}, expected ${b}`); };
  near(K.wStar(1, 4, 0.25), 0.875, "PP4 w*");
  near(K.varStar(1, 4, 0.25), 0.9375, "PP4 Var at w*");
  near(K.varC(0.875, 1, 4, 0.25), 0.9375, "PP4 Var(e_c)(w*)");
  near(K.varC(0.5, 1, 4, 0.25), 1.5, "PP4 simple average");
  near(K.wInv(1, 4), 0.8, "PP4 inverse-variance w");
  near(K.varC(0.8, 1, 4, 0.25), 0.96, "PP4 Var at inverse-variance w");
  near(K.rhoThreshold(1, 4), 0.5, "threshold rho_{1,2} at ratio 4");
  { const b = K.beatInterval(1, 4, 0.5); near(b[0], 1, "beat interval lo at threshold"); near(b[1], 1, "beat interval hi at threshold"); }
  near(K.wStar(1, 1.5, 0), K.wInv(1, 1.5), "zero covariance: w* = inverse-variance w");
  const checkFailed = problems.length > 0;
  if (checkFailed) console.error(`makeCombo2[${uid}] self-check FAILED:\n  ` + problems.join("\n  "));

  // ---- formatting helpers -------------------------------------------------
  const MINUS = "−";
  const fmt = (x, d) => {
    let s = Number(x).toFixed(d);
    if (/^-0(\.0*)?$/.test(s)) s = s.slice(1);
    return s.replace("-", MINUS);
  };
  const pct = x => (x > 0 ? "+" : "") + fmt(100 * x, 1) + "%";
  const offGrid = x => Math.abs(x * 100 - Math.round(x * 100)) > 1e-6;
  const esc = s => String(s).replace(/[&<>"]/g, ch => ({ "&": "&amp;", "<": "&lt;", ">": "&gt;", '"': "&quot;" }[ch]));
  // MathML pieces (course notation)
  const eSub = s => `<msub><mi>e</mi><mrow>${s}</mrow></msub>`;
  const e1 = eSub(`<mn>1</mn><mo>,</mo><mi>t</mi>`);
  const e2 = eSub(`<mn>2</mn><mo>,</mo><mi>t</mi>`);
  const ec = eSub(`<mi>c</mi><mo>,</mo><mi>t</mi>`);
  const VarM = x => `<mtext>Var</mtext><mo stretchy="false">(</mo>${x}<mo stretchy="false">)</mo>`;
  const CovM = `<mtext>Cov</mtext><mo stretchy="false">(</mo>${e1}<mo>,</mo>${e2}<mo stretchy="false">)</mo>`;
  const m = s => `<math>${s}</math>`;
  // rho_{1,2}: comma-separated two-index subscript (a bare 12 reads as the lag-12 autocorrelation)
  const rho12 = `<msub><mi>ρ</mi><mrow><mn>1</mn><mo>,</mo><mn>2</mn></mrow></msub>`;
  const wS = `<msup><mi>w</mi><mo>*</mo></msup>`;
  const wI = `<mi>w</mi>`;
  const V1m = m(VarM(e1)), V2m = m(VarM(e2)), VCm = m(VarM(ec)), COVm = m(CovM), RHOm = m(rho12), WSm = m(wS), Wm = m(wI);

  // ---- current derived quantities (shared by render and the bubbles) ------------
  const derive = () => {
    const v1 = V1, v2 = ratio * V1;
    return { v1, v2, ws: K.wStar(v1, v2, rho), vs: K.varStar(v1, v2, rho), vw: K.varC(w, v1, v2, rho),
             wi: K.wInv(v1, v2), thr: K.rhoThreshold(v1, v2), cov: K.cov(v1, v2, rho) };
  };

  // ---- bubble texts: one or two sentences in course notation, plus a ref key --------
  const TIPS = {
    w: { ref: "w", html: d => `${Wm} is the weight on forecast 1; forecast 2 gets 1 − ${Wm}, so the weights sum to one. ` +
         `Unhatted ${Wm} is a generic, fixed weight; below 0 or above 1 one of the forecasts is shorted.` },
    rho12: { ref: "rho12", html: d => `${RHOm} is the correlation between the two forecasters' errors, ${COVm} / √(${V1m} ${V2m}). ` +
         `The comma-separated pair of forecaster indices keeps it apart from the lag-<i>k</i> autocorrelation ${m("<msub><mi>ρ</mi><mi>k</mi></msub>")}.` },
    var_ratio: { ref: "var_ratio", html: d => `How much worse forecast 2 is: ${V2m} / ${V1m}, with ${V1m} held at 1. ` +
         `At 1 or above, forecast 1 is always the better forecast.` },
    to_wstar: { ref: "wstar", html: d => (wsOff(d.ws)
         ? `${WSm} = ${fmt(d.ws, 3)} lies beyond the ${Wm} axis [${fmt(WMIN, 2)}, ${fmt(WMAX, 2)}]: the button <b>clamps</b> ${Wm} to ${fmt(clampW(d.ws), 2)}, not ${WSm}. `
         : `Moves ${Wm} to ${WSm} = ${fmt(d.ws, 3)}, the weight that minimizes ${VCm} (Enders 2.73). `) +
         `Starred because it uses the true variances and covariance; weights estimated from data carry hats.` },
    to_threshold: { ref: "threshold", html: d => thrAvail()
         ? `Sets ${RHOm} to √(${V1m} / ${V2m}) = ${fmt(d.thr, 3)}, where ${COVm} = ${V1m}. There ${WSm} = 1 and no weight in [0, 1] beats forecast 1.`
         : ratio > 1
         ? `Unavailable at this variance ratio: the threshold √(${V1m} / ${V2m}) = ${fmt(d.thr, 3)} is above the slider's ${fmt(G.rho.max, 2)}. Every ${RHOm} on the slider is below it, so ${WSm} &lt; 1 here.`
         : `Unavailable at equal variances: the threshold would be ${RHOm} = 1, outside the slider. With ${V1m} = ${V2m}, every ${RHOm} &lt; 1 gives ${WSm} = ½.` },
    preset_uncorrelated: { ref: "inv_var", html: d => `Sets ${RHOm} = 0 at variance ratio 1.5. With zero covariance ${WSm} equals the inverse-variance weight: the worse forecast still gets weight, just less.` },
    preset_fig09: { ref: "short", html: d => `The notes' §7.5 figure at ${RHOm} = 0.9 (variance ratio 1.5). Past the threshold ${WSm} &gt; 1: the optimum shorts forecast 2.` },
    preset_negative: { ref: "rho12", html: d => `Sets ${RHOm} = −0.5 at variance ratio 1.5: errors that tend to offset, so combining gains the most.` },
    reset: { ref: "varec", html: d => `Restores the starting ${Wm} = ${fmt(DEF.w, 2)}, ${RHOm} = ${fmt(DEF.rho, 2)} and variance ratio ${fmt(DEF.ratio, 2)}.` },
    f1: { ref: "varec", html: d => `Forecast 1 used alone, ${Wm} = 1: ${VCm} = ${V1m} = 1. It is the better of the two forecasts.` },
    f2: { ref: "var_ratio", html: d => `Forecast 2 used alone, ${Wm} = 0: ${VCm} = ${V2m} = ${fmt(d.v2, 2)}.` },
    avg: { ref: "simple_avg", html: d => `The simple average, ${Wm} = ½, which is fixed, not estimated, so it carries no hat. Here ${VCm} = ${fmt(K.varC(0.5, d.v1, d.v2, rho), 3)}; it beats forecast 1 only when its diamond is below the dashed line.` },
    inv: { ref: "inv_var", html: d => `Inverse-variance weight ${V1m}⁻¹ / (${V1m}⁻¹ + ${V2m}⁻¹) = ${fmt(d.wi, 3)}: the ${WSm} you would get if ${COVm} were 0. It ignores ${RHOm}; §7.6 Method 2 estimates it from MSEs.` },
    wstar: { ref: "wstar", html: d => `${WSm} = (${V2m} − ${COVm}) / (${V1m} + ${V2m} − 2${COVm}) = ${fmt(d.ws, 3)}, the bottom of the curve (Enders 2.73, Bates-Granger). ` +
         `${VCm} there is ${fmt(d.vs, 3)}.` },
    wcur: { ref: "w", html: d => `Your weight, ${Wm} = ${fmt(w, offGrid(w) ? 3 : 2)}, with ${VCm} = ${fmt(d.vw, 3)}. Drag the ${Wm} slider to move it along the curve.` },
    beat: { ref: "threshold", html: d => `Weights at which the combination beats the better forecast alone: ${VCm} &lt; ${V1m}. ` +
         `Inside [0, 1] this region exists only while ${COVm} &lt; ${V1m}.` },
    refline: { ref: "varec", html: d => `${V1m} = 1, the better forecast used on its own. Any combination whose curve dips below this line beats it.` },
    r_varw: { ref: "varec", html: d => `${VCm} = ${Wm}²${V1m} + (1 − ${Wm})²${V2m} + 2${Wm}(1 − ${Wm})${COVm}, evaluated at your ${Wm}; the percent compares it with the better forecast, ${V1m} = 1.` },
    r_wstar: { ref: "wstar", html: d => `The optimal population weight on forecast 1 (Enders 2.73). It can leave [0, 1]: above 1 it shorts forecast 2.` },
    r_varstar: { ref: "wstar", html: d => `${VCm} at ${WSm}, the smallest combined error variance any weight can reach: (${V1m}${V2m} − ${COVm}²) / (${V1m} + ${V2m} − 2${COVm}).` },
    r_threshold: { ref: "threshold", html: d => `The ${RHOm} at which ${COVm} = ${V1m}, namely √(${V1m} / ${V2m}). Above it, ${WSm} &gt; 1 and only shorting beats forecast 1.` +
         (ratio > 1 && !thrAvail() ? ` At this ratio it is ${fmt(d.thr, 3)}, above the slider's ${fmt(G.rho.max, 2)}.` : "") }
  };

  // ---- DOM kit ------------------------------------------------------------
  const el = (tag, cls, html) => { const n = document.createElement(tag); if (cls) n.className = cls; if (html != null) n.innerHTML = html; return n; };
  const button = (label, action, cls, tip) => { const b = document.createElement("button"); b.type = "button"; b.className = "cx2-btn" + (cls ? " " + cls : ""); b.innerHTML = label; b.dataset.action = action; if (tip) b.dataset.tip = tip; return b; };
  const I_MARK = `<span class="cx2-i" aria-hidden="true">ⓘ</span>`;
  const root = el("div", "cx2");
  root.dataset.widget = "combo2";
  root.dataset.uid = String(uid);

  const mkSlider = (labelHTML, g, val, action, aria, tip) => {
    const ctl = el("label", "cx2-ctl");
    ctl.dataset.tip = tip;
    const s = document.createElement("input");
    s.type = "range"; s.min = g.min; s.max = g.max; s.step = g.step; s.value = val; s.dataset.action = action;
    s.setAttribute("aria-label", aria);
    const out = el("span", "cx2-cval");
    ctl.append(el("span", "cx2-clab", labelHTML), s, out);
    return { ctl, s, out };
  };
  const wC = mkSlider(m(wI) + " =", G.w, w, "w", "weight w on forecast 1", "w");
  const rC = mkSlider(m(rho12) + " =", G.rho, rho, "rho", "error correlation rho sub 1, 2", "rho12");
  const vC = mkSlider(m(`${VarM(e2)}<mo>/</mo>${VarM(e1)}`) + " =", G.ratio, ratio, "ratio",
                      "variance ratio, Var of e 2 over Var of e 1", "var_ratio");
  const bar1 = el("div", "cx2-bar");
  bar1.append(wC.ctl, rC.ctl, vC.ctl);

  const bar2 = el("div", "cx2-bar");
  const toStarB = button(`${m(wI)} → ${m(wS)}`, "w-to-wstar", null, "to_wstar");
  const toThrB = button(`${m(rho12)} → threshold`, "rho-to-threshold", null, "to_threshold");
  const grpA = el("div", "cx2-grp"); grpA.append(el("span", "cx2-lab", "snap"), toStarB, toThrB);
  bar2.appendChild(grpA);
  const PRESETS = [
    { label: `uncorrelated`, action: "preset-uncorrelated", tip: "preset_uncorrelated", set: { rho: 0, ratio: 1.5 } },
    { label: `§7.5 figure, ${m(rho12)} = 0.9`, action: "preset-fig09", tip: "preset_fig09", set: { rho: 0.9, ratio: 1.5 } },
    { label: `negative ${m(rho12)}`, action: "preset-negative", tip: "preset_negative", set: { rho: -0.5, ratio: 1.5 } }
  ];
  const presetBtns = [];
  if (SHOW_PRESETS) {
    const grpP = el("div", "cx2-grp"); grpP.appendChild(el("span", "cx2-lab", "try"));
    PRESETS.forEach(p => { const b = button(p.label, p.action, null, p.tip); presetBtns.push([b, p]); grpP.appendChild(b); });
    bar2.appendChild(grpP);
  }
  const resetB = button("Reset", "reset", "cx2-ghost", "reset");
  resetB.style.marginLeft = "auto";
  bar2.appendChild(resetB);
  if (checkFailed) bar2.appendChild(el("span", "cx2-warnline", "⚠ self-check failed — see console"));

  const msg = el("div", "cx2-msg");
  // readout: static structure (so a bubble anchored on a field survives a re-render); values
  // are written into the data-v spans
  const read = el("div", "cx2-read");
  const fld = (tip, html) => `<span class="cx2-f" data-tip="${tip}" tabindex="0">${html}${I_MARK}</span>`;
  read.innerHTML =
    fld("r_varw", `${VCm} at ${Wm} = <span data-v="w"></span>: <b data-v="vw"></b> <span class="cx2-dim">(<span data-v="vwp"></span> vs better)</span>`) +
    `<span class="cx2-dim"> · </span>` + fld("r_wstar", `${WSm} = <b data-v="ws"></b>`) +
    `<span class="cx2-dim"> · </span>` + fld("r_varstar", `at ${WSm}: <b data-v="vs"></b> <span class="cx2-dim">(<span data-v="vsp"></span> vs better)</span>`) +
    `<span class="cx2-dim"> · </span>` + fld("r_threshold", `<span class="cx2-dim">threshold ${RHOm} = <span data-v="thr"></span></span>`);
  const RV = {}; read.querySelectorAll("[data-v]").forEach(n => { RV[n.dataset.v] = n; });
  const chartBox = el("div", "cx2-chart");
  const legend = el("div", "cx2-legend");
  root.append(bar1, bar2, msg, read, chartBox, legend);

  // ---- geometry -------------------------------------------------------------
  // Deck: viewBox 1800 x 570, mounted at about 1800 CSS px on a 1920-wide slide, so one unit
  // is about one screen px and FS = 26 matches the HTML type floor. Notes: viewBox 1000 x 460
  // in a ~800 px column (0.8 px per unit), FS = 24 renders at about 19 px.
  const W = NOTES ? 1000 : 1800, H = NOTES ? 460 : 570, FS = NOTES ? 24 : 26;
  const L = NOTES ? 100 : 130, R = W - (NOTES ? 22 : 30), PT = 46, PB = H - 70;
  const WMIN = G.w.min, WMAX = G.w.max;
  const sx = x => L + (x - WMIN) / (WMAX - WMIN) * (R - L);
  const col = { ink: "#16202b", muted: "#5d6b78", line: "#d7dee5", soft: "#e7ecf1", zero: "#aab4bf",
                grey: "#8a94a0", blue: "#2c6fbb", orange: "#e4572e", navy: "#0b1f3a", green: "#2f8f6b",
                purple: "#6b4fa0", greenFill: "rgba(47,143,107,0.22)", shortFill: "#f4f6f8" };
  const P = (x, y) => `${x.toFixed(1)},${y.toFixed(1)}`;
  // Vertical scale: discrete tops so it does not breathe while a slider moves, capped at 2.5 so
  // the region near the minimum (all the action is at or below Var(e_1) = 1) stays legible.
  // Forecast 2 (Var(e_2) up to 6) and the curve's ends go off the top; markers become arrows.
  const TOPS = [1.6, 2, 2.5];
  const yTop = r => { const need = Math.max(1.25, 1.12 * Math.min(r, 2.2)); for (const t of TOPS) if (t >= need) return t; return TOPS[TOPS.length - 1]; };
  const yStep = t => t <= 2 ? 0.25 : 0.5;
  // SVG text with a subscript, e.g. Var(e_{c,t}): tspans with dy (portable, no baseline-shift)
  const subT = (base, s, after) => `${base}<tspan dy="${(FS * 0.3).toFixed(1)}" font-size="${(FS * 0.75).toFixed(1)}">${s}</tspan><tspan dy="${(-FS * 0.3).toFixed(1)}">${after}</tspan>`;

  const clipId = `cx2clip${uid}`;
  let xt = "";
  const NTICK = Math.round((WMAX - WMIN) / 0.25);        // 0.25 spacing over the whole w axis
  for (let k = 0; k <= NTICK; k++) {
    const v = WMIN + k * 0.25, X = sx(v);
    const isInt = Math.abs(v - Math.round(v)) < 1e-9, major = Math.abs(v * 2 - Math.round(v * 2)) < 1e-9;
    const digits = isInt ? 0 : (major ? 1 : 2);
    xt += `<line x1="${X.toFixed(1)}" y1="${PB}" x2="${X.toFixed(1)}" y2="${PB + 8}" stroke="${col.grey}"/>`;
    if (major || !NOTES) xt += `<text x="${X.toFixed(1)}" y="${PB + 34}" text-anchor="middle" font-size="${FS}" fill="${col.muted}">${fmt(v, digits)}</text>`;
  }
  chartBox.innerHTML =
`<svg viewBox="0 0 ${W} ${H}" role="img" aria-label="Variance of the combined forecast error as a function of the weight w on forecast 1, with the better forecast's variance as a dashed reference line">
  <defs><clipPath id="${clipId}"><rect x="${L}" y="${PT}" width="${R - L}" height="${PB - PT}"/></clipPath></defs>
  <rect x="${sx(WMIN).toFixed(1)}" y="${PT}" width="${(sx(0) - sx(WMIN)).toFixed(1)}" height="${PB - PT}" fill="${col.shortFill}"/>
  <rect x="${sx(1).toFixed(1)}" y="${PT}" width="${(sx(WMAX) - sx(1)).toFixed(1)}" height="${PB - PT}" fill="${col.shortFill}"/>
  <text x="${((sx(WMIN) + sx(0)) / 2).toFixed(1)}" y="${PT - 14}" text-anchor="middle" font-size="${FS}" fill="${col.muted}">${NOTES ? "shorts forecast 1" : "<tspan font-style=\"italic\">w</tspan> &lt; 0: shorts forecast 1"}</text>
  <text x="${((sx(1) + sx(WMAX)) / 2).toFixed(1)}" y="${PT - 14}" text-anchor="middle" font-size="${FS}" fill="${col.muted}">${NOTES ? "shorts forecast 2" : "<tspan font-style=\"italic\">w</tspan> &gt; 1: shorts forecast 2"}</text>
  <g class="cx2-gy"></g>
  ${xt}
  <text x="${((L + R) / 2).toFixed(1)}" y="${H - 6}" text-anchor="middle" font-size="${FS}" fill="${col.muted}"><tspan font-style="italic">w</tspan> = weight on forecast 1 (forecast 2 gets 1 − <tspan font-style="italic">w</tspan>)</text>
  <text transform="translate(${(FS * 1.1).toFixed(1)} ${((PT + PB) / 2).toFixed(1)}) rotate(-90)" text-anchor="middle" font-size="${FS}" font-weight="600" fill="${col.ink}">${subT("Var(e", "c,t", ")")}</text>
  <g class="cx2-gdyn" clip-path="url(#${clipId})"></g>
  <g class="cx2-gtop"></g>
  <line x1="${L}" y1="${PB}" x2="${R}" y2="${PB}" stroke="${col.ink}" stroke-width="1.6"/>
  <line x1="${L}" y1="${PT}" x2="${L}" y2="${PB}" stroke="${col.ink}" stroke-width="1.2"/>
</svg>`;
  const svg = chartBox.querySelector("svg");
  const gY = svg.querySelector(".cx2-gy"), gDyn = svg.querySelector(".cx2-gdyn"), gTop = svg.querySelector(".cx2-gtop");

  // legend (HTML, so it reads at the type floor; each entry is the keyboard route to its glyph's bubble)
  const sw = inner => `<svg class="cx2-sw" viewBox="0 0 40 24" aria-hidden="true">${inner}</svg>`;
  const star = (cx, cy, r, fill) => {
    let pts = [];
    for (let i = 0; i < 10; i++) { const a = -Math.PI / 2 + i * Math.PI / 5, rr = i % 2 ? r * 0.45 : r; pts.push(P(cx + rr * Math.cos(a), cy + rr * Math.sin(a))); }
    return `<polygon points="${pts.join(" ")}" fill="${fill}" stroke="#fff" stroke-width="1.5"/>`;
  };
  const diamond = (cx, cy, r, fill) => `<polygon points="${P(cx, cy - r)} ${P(cx + r, cy)} ${P(cx, cy + r)} ${P(cx - r, cy)}" fill="${fill}" stroke="#fff" stroke-width="1.5"/>`;
  const square = (cx, cy, r, fill) => `<rect x="${(cx - r).toFixed(1)}" y="${(cy - r).toFixed(1)}" width="${(2 * r).toFixed(1)}" height="${(2 * r).toFixed(1)}" fill="${fill}" stroke="#fff" stroke-width="1.5"/>`;
  const li = (tip, swatch, html) => `<span class="cx2-li" data-tip="${tip}" tabindex="0">${sw(swatch)}<span>${html}</span>${I_MARK}</span>`;
  legend.innerHTML =
    li("f1", `<circle cx="20" cy="12" r="8" fill="${col.navy}"/>`, `forecast 1 alone (<i>w</i> = 1)`) +
    li("f2", `<circle cx="20" cy="12" r="8" fill="${col.orange}"/>`, `forecast 2 alone (<i>w</i> = 0)`) +
    li("avg", diamond(20, 12, 9, col.grey), `simple average (<i>w</i> = ½)`) +
    li("inv", square(20, 12, 7, col.purple), `inverse-variance`) +
    li("wstar", star(20, 12, 11, col.green), `${WSm} (Bates-Granger)`) +
    li("wcur", `<line x1="20" y1="0" x2="20" y2="24" stroke="${col.blue}" stroke-width="3"/><circle cx="20" cy="12" r="6" fill="${col.blue}"/>`, `your ${Wm}`) +
    li("refline", `<line x1="2" y1="12" x2="38" y2="12" stroke="${col.navy}" stroke-width="2.2" stroke-dasharray="7 5"/>`, `better forecast`) +
    li("beat", `<rect x="2" y="4" width="36" height="16" fill="${col.greenFill}" stroke="${col.green}"/>`, `beats it`);

  // ---- render ----------------------------------------------------------------
  const render = () => {
    const d = derive();
    const { v1, v2, ws, vs, vw, wi, thr } = d;
    const top = yTop(ratio), step = yStep(top);
    const sy = v => PT + (top - v) / top * (PB - PT);
    const best = Math.min(v1, v2);
    const atThr = ratio > 1 && Math.abs(rho - thr) <= AT_TOL;
    const regime = atThr ? "threshold" : (ws > 1 ? "short" : "pays");

    // slider readouts and the threshold button's state
    wC.s.value = w; rC.s.value = rho; vC.s.value = ratio;
    wC.out.textContent = fmt(w, offGrid(w) ? 3 : 2);
    rC.out.textContent = fmt(rho, offGrid(rho) ? 3 : 2);
    vC.out.textContent = fmt(ratio, 2);
    toThrB.setAttribute("aria-disabled", thrAvail() ? "false" : "true");

    // y grid
    let yg = "";
    for (let v = 0; v <= top + 1e-9; v += step) {
      const Y = sy(v);
      yg += `<line x1="${L}" y1="${Y.toFixed(1)}" x2="${R}" y2="${Y.toFixed(1)}" stroke="${col.soft}"/>` +
            `<text x="${L - 10}" y="${(Y + FS * 0.35).toFixed(1)}" text-anchor="end" font-size="${FS}" fill="${col.muted}">${fmt(v, step < 0.5 ? 2 : 1)}</text>`;
    }
    yg += `<line x1="${sx(0).toFixed(1)}" y1="${PT}" x2="${sx(0).toFixed(1)}" y2="${PB}" stroke="${col.zero}" stroke-width="1.2"/>` +
          `<line x1="${sx(1).toFixed(1)}" y1="${PT}" x2="${sx(1).toFixed(1)}" y2="${PB}" stroke="${col.zero}" stroke-width="1.2"/>`;
    gY.innerHTML = yg;

    // curve, sampled densely (closed form; sampling is only for drawing)
    // The geometry itself is cut at the plot top (not only clipped visually): with the wide w
    // axis the curve's ends rise far above Var = top, and an unclipped polyline's bounding box
    // would reach above the slide. Pieces above the top are dropped; each crossing gets an
    // exact end point on y = PT (linear interpolation between samples), as separate subpaths.
    const N = 400, segs = [];
    let cur = null, xp = null, vp = null;
    for (let i = 0; i <= N; i++) {
      const x = WMIN + (WMAX - WMIN) * i / N, v = K.varC(x, v1, v2, rho), inside = v <= top;
      if (i > 0 && (vp <= top) !== inside) {                    // crossing between samples
        const xc = xp + (x - xp) * (top - vp) / (v - vp);
        if (inside) cur = [P(sx(xc), PT)]; else { cur.push(P(sx(xc), PT)); segs.push(cur); cur = null; }
      }
      if (inside) { if (!cur) cur = []; cur.push(P(sx(x), sy(v))); }
      xp = x; vp = v;
    }
    if (cur && cur.length) segs.push(cur);
    const curveD = segs.filter(s => s.length > 1).map(s => "M" + s.join(" L")).join(" ");
    // beat region: analytic interval between the roots of Var(e_c)(w) = Var(e_1)
    let fill = "";
    const bi = K.beatInterval(v1, v2, rho);
    const lo = Math.max(WMIN, bi[0]), hi = Math.min(WMAX, bi[1]);
    if (!atThr && hi - lo > 1e-6) {
      const fp = [], M = 160;
      for (let i = 0; i <= M; i++) { const x = lo + (hi - lo) * i / M; fp.push(P(sx(x), sy(K.varC(x, v1, v2, rho)))); }
      fp.push(P(sx(hi), sy(best)), P(sx(lo), sy(best)));
      fill = `<polygon class="cx2-beat" data-tip="beat" points="${fp.join(" ")}" fill="${col.greenFill}" stroke="none"/>`;
    }
    const refY = sy(best);
    gDyn.innerHTML =
      fill +
      `<line x1="${L}" y1="${refY.toFixed(1)}" x2="${R}" y2="${refY.toFixed(1)}" stroke="${col.navy}" stroke-width="2.2" stroke-dasharray="12 8"/>` +
      `<line data-tip="refline" x1="${L}" y1="${refY.toFixed(1)}" x2="${R}" y2="${refY.toFixed(1)}" stroke="transparent" stroke-width="16"/>` +
      `<path class="cx2-curve" d="${curveD}" fill="none" stroke="${col.ink}" stroke-width="3.5" stroke-linejoin="round" pointer-events="none"/>`;

    // markers and labels (unclipped layer). A marker above the top becomes an up-arrow at the
    // top edge carrying its value.
    const mk = [];
    const refLab = NOTES ? `forecast 1 alone: ${subT("Var(e", "1,t", ")")} = 1` : `the better forecast alone: ${subT("Var(e", "1,t", ")")} = 1`;
    // left of w = 0 the curve is always above Var(e_1) (w* > 0 whenever Var(e_2) >= Var(e_1)),
    // so the space just under the dashed line there is free for its label
    mk.push(`<text data-tip="refline" x="${(L + 10).toFixed(1)}" y="${(refY + FS + 6).toFixed(1)}" font-size="${FS}" fill="${col.navy}" paint-order="stroke" stroke="#fff" stroke-width="5">${refLab}</text>`);
    const xw = sx(w), yw = sy(vw);
    mk.push(`<line class="cx2-wline" data-tip="wcur" x1="${xw.toFixed(1)}" y1="${PT}" x2="${xw.toFixed(1)}" y2="${PB}" stroke="${col.blue}" stroke-width="2.5" opacity="0.85"/>`);
    const upArrow = (x, color, tip, text, rowOffset) => {
      const anc = x > R - 220 ? "end" : (x < L + 220 ? "start" : "middle");
      const dx = anc === "end" ? -16 : (anc === "start" ? 16 : 0);
      const ty = PT + 22 + FS + (rowOffset || 0);
      return `<g data-tip="${tip}" class="cx2-arrow"><polygon points="${P(x, PT + 2)} ${P(x - 11, PT + 20)} ${P(x + 11, PT + 20)}" fill="${color}"/>` +
             `<text x="${(x + dx).toFixed(1)}" y="${ty.toFixed(1)}" text-anchor="${anc}" font-size="${FS}" font-weight="700" fill="${color}" paint-order="stroke" stroke="#fff" stroke-width="6">${text}</text></g>`;
    };
    const glyph = (x, v, tip, draw, color, arrowText, row) =>
      sy(v) >= PT ? `<g data-tip="${tip}">${draw(x, sy(v))}</g>` : upArrow(x, color, tip, arrowText, row);
    const vAvg = K.varC(0.5, v1, v2, rho), vInv = K.varC(wi, v1, v2, rho);
    mk.push(glyph(sx(0.5), vAvg, "avg", (x, y) => diamond(x, y, 12, col.grey), col.grey, `Var at average = ${fmt(vAvg, 2)}`, FS + 6));
    mk.push(glyph(sx(wi), vInv, "inv", (x, y) => square(x, y, 9, col.purple), col.purple, `Var at inverse-variance = ${fmt(vInv, 2)}`, 2 * (FS + 6)));
    mk.push(glyph(sx(0), v2, "f2", (x, y) => `<circle cx="${x.toFixed(1)}" cy="${y.toFixed(1)}" r="10" fill="${col.orange}" stroke="#fff" stroke-width="2"/>`, col.orange, `Var of forecast 2 = ${fmt(v2, 2)}`));
    mk.push(`<g data-tip="f1"><circle cx="${sx(1).toFixed(1)}" cy="${sy(v1).toFixed(1)}" r="10" fill="${col.navy}" stroke="#fff" stroke-width="2"/></g>`);
    // w*
    if (ws <= WMAX) {
      const xs = sx(ws), ys = sy(vs);
      const anchor = xs > R - 120 ? "end" : (xs < L + 120 ? "start" : "middle");
      // label on the axis below the star when there is room; above the star when Var(w*) is
      // near zero (the star sits on the axis and would collide with the label)
      const below = ys < PB - 2.2 * FS;
      const ly = below ? PB - 10 : ys - 15 - 0.5 * FS;
      mk.push(`<line x1="${xs.toFixed(1)}" y1="${ys.toFixed(1)}" x2="${xs.toFixed(1)}" y2="${PB}" stroke="${col.green}" stroke-width="2" stroke-dasharray="4 5" pointer-events="none"/>`);
      mk.push(`<g data-tip="wstar">${star(xs, ys, 15, col.green)}` +
              `<text class="cx2-wslab" x="${xs.toFixed(1)}" y="${ly.toFixed(1)}" text-anchor="${anchor}" font-size="${FS}" font-weight="700" fill="${col.green}" paint-order="stroke" stroke="#fff" stroke-width="6">w* = ${fmt(ws, 3)}</text></g>`);
    } else {
      // fallback only: the controls keep w* <= 2.101 < WMAX (see G.w); reachable via the QA hook
      const ya = PT + 40;                    // the curve is low at the right edge when w* > WMAX
      mk.push(`<g data-tip="wstar"><polygon points="${P(R - 2, ya)} ${P(R - 22, ya - 12)} ${P(R - 22, ya + 12)}" fill="${col.green}"/>` +
              `<text x="${(R - 28).toFixed(1)}" y="${(ya + FS * 0.35).toFixed(1)}" text-anchor="end" font-size="${FS}" font-weight="700" fill="${col.green}" paint-order="stroke" stroke="#fff" stroke-width="6">w* = ${fmt(ws, 3)} (off chart)</text></g>`);
    }
    mk.push(glyph(xw, vw, "wcur", (x, y) => `<circle class="cx2-wdot" cx="${x.toFixed(1)}" cy="${y.toFixed(1)}" r="9" fill="${col.blue}" stroke="#fff" stroke-width="2"/>`, col.blue, `Var at your w = ${fmt(vw, 2)}`, 3 * (FS + 6)));
    gTop.innerHTML = mk.join("");

    // regime message
    const gain = vs / best - 1;
    if (regime === "pays") {
      msg.innerHTML = `<span class="cx2-tag cx2-good">diversification pays</span>` +
        `<span>${WSm} = ${fmt(ws, 3)} mixes both forecasts; the best mix cuts ${VCm} ${fmt(-100 * gain, 1)}% below the better forecast alone. ` +
        `<span class="cx2-dim">(${COVm} = ${fmt(d.cov, 3)} &lt; ${V1m} = 1)</span></span>`;
    } else if (regime === "threshold") {
      msg.innerHTML = `<span class="cx2-tag cx2-flat">no gain</span>` +
        `<span>${RHOm} at threshold ${m(`<msqrt><mn>1</mn><mo>/</mo><mn>${fmt(ratio, 2)}</mn></msqrt>`)} = ${fmt(thr, 3)}: ${COVm} = ${V1m}, ${WSm} = 1, ` +
        `so forecast 1 alone is the best combination.</span>`;
    } else {
      msg.innerHTML = `<span class="cx2-tag cx2-short">${WSm} &gt; 1</span>` +
        `<span>optimal combination shorts the worse forecast: ${WSm} = ${fmt(ws, 3)} puts ${fmt(1 - ws, 3)} on forecast 2. ` +
        `<span class="cx2-dim">No weight in [0, 1] beats forecast 1 (${COVm} = ${fmt(d.cov, 3)} &gt; ${V1m}).</span></span>`;
    }
    root.dataset.regime = regime;

    // readout values
    RV.w.textContent = fmt(w, offGrid(w) ? 3 : 2);
    RV.vw.textContent = fmt(vw, 3); RV.vwp.textContent = pct(vw / best - 1);
    RV.ws.textContent = fmt(ws, 3) + (wsOff(ws) ? ` (beyond the axis; w → w* clamps to ${fmt(clampW(ws), 2)})` : "");
    RV.vs.textContent = fmt(vs, 3); RV.vsp.textContent = pct(gain);
    RV.thr.textContent = ratio > 1 ? fmt(thr, 3) + (thrAvail() ? "" : ` (above the slider's ${fmt(G.rho.max, 2)})`) : "none (equal variances)";

    refreshTip();
  };

  // ---- info bubbles ----------------------------------------------------------------
  // One bubble per instance. It is moved in the DOM to sit right after its trigger (after the
  // whole chart for SVG glyphs), so Tab goes from a trigger into the bubble's link; it is
  // positioned absolutely against the root, converting screen px to local px through the
  // reveal.js scale.
  const tip = el("div", "cx2-tip");
  tip.id = `cx2tip${uid}`; tip.setAttribute("role", "tooltip"); tip.hidden = true;
  root.appendChild(tip);
  let openTrig = null, openKey = null, pinned = false, hideTimer = null, quiet = false;
  const isControl = t => t.tagName === "BUTTON" || t.tagName === "LABEL";
  const tipHTML = key => {
    const t = TIPS[key]; if (!t) return "";
    const r = REFS[t.ref] || { label: "Notes §7.5" };
    const lab = esc(r.label || "");
    const ptr = r.href ? `<a href="${esc(r.href)}">${lab}</a>` : lab;
    return `<div class="cx2-tip-b">${t.html(derive())}</div><div class="cx2-tip-ref">Defined on: ${ptr}</div>`;
  };
  const place = trig => {
    const rr = root.getBoundingClientRect();
    const k = rr.width / (root.offsetWidth || 1) || 1;
    const tr = trig.getBoundingClientRect();
    tip.style.visibility = "hidden"; tip.hidden = false; tip.style.left = "0px"; tip.style.top = "0px";
    const tw = tip.offsetWidth, th = tip.offsetHeight, RW = root.offsetWidth, RH = root.offsetHeight;
    let x = (tr.left - rr.left) / k;
    x = Math.max(0, Math.min(RW - tw, x));
    const tTop = (tr.top - rr.top) / k, tBot = (tr.bottom - rr.top) / k;
    let y = tBot + 8;
    if (y + th > RH && tTop - th - 8 >= 0) y = tTop - th - 8;
    tip.style.left = x.toFixed(1) + "px"; tip.style.top = y.toFixed(1) + "px";
    tip.style.visibility = "";
  };
  const show = (trig, pin) => {
    clearTimeout(hideTimer);
    const key = trig.dataset.tip;
    if (!TIPS[key]) return;
    if (openTrig && openTrig !== trig) openTrig.removeAttribute("aria-describedby");
    openTrig = trig; openKey = key; pinned = !!pin;
    tip.innerHTML = tipHTML(key);
    tip.dataset.key = key;
    const anchor = trig.closest("svg") ? chartBox : trig;
    if (anchor.nextSibling !== tip) anchor.after(tip);
    trig.setAttribute("aria-describedby", tip.id);
    place(trig);
  };
  const hide = () => {
    clearTimeout(hideTimer);
    if (openTrig) openTrig.removeAttribute("aria-describedby");
    tip.hidden = true; openTrig = null; openKey = null; pinned = false;
  };
  const scheduleHide = () => { clearTimeout(hideTimer); hideTimer = setTimeout(() => { if (!pinned) hide(); }, 220); };
  // after a render: SVG glyphs are rebuilt, so re-find the trigger by key; refresh text (it
  // quotes live values) and position
  function refreshTip() {
    if (!openTrig || tip.hidden) return;
    if (!root.contains(openTrig)) {
      const again = root.querySelector(`.cx2-chart [data-tip="${openKey}"]`);
      if (!again) { hide(); return; }
      openTrig = again;
    }
    tip.innerHTML = tipHTML(openKey);
    place(openTrig);
  }
  const trigOf = e => { const t = e.target.closest && e.target.closest("[data-tip]"); return t && root.contains(t) && !tip.contains(t) ? t : null; };
  root.addEventListener("pointerover", e => {
    if (tip.contains(e.target)) { clearTimeout(hideTimer); return; }
    const t = trigOf(e);
    if (!t || e.pointerType === "touch") return;
    if (t !== openTrig) { if (!pinned) show(t, false); } else clearTimeout(hideTimer);
  });
  root.addEventListener("pointerout", e => {
    if (!openTrig || pinned) return;
    const to = e.relatedTarget;
    if (to && (tip.contains(to) || openTrig.contains(to))) return;
    scheduleHide();
  });
  // touch: tap opens (controls still act); tap on a non-control trigger toggles
  let lastDown = { type: null, t: 0 };
  root.addEventListener("pointerdown", e => {
    lastDown = { type: e.pointerType, t: Date.now() };
    if (tip.contains(e.target)) return;
    const t = trigOf(e);
    if (e.pointerType === "touch") {
      if (t) { if (t === openTrig && !isControl(t)) hide(); else show(t, true); }
      else if (openTrig) hide();
      return;
    }
    // mouse press on a slider: hide so the bubble does not cover the chart while dragging
    if (t && t.tagName === "LABEL" && t === openTrig && !pinned) hide();
  });
  // click (mouse) on a non-control trigger pins the bubble so its link can be reached
  root.addEventListener("click", e => {
    if (e.pointerType === "touch" || (lastDown.type === "touch" && Date.now() - lastDown.t < 800)) return;
    const t = trigOf(e);
    if (!t || isControl(t)) return;
    if (t === openTrig && pinned) hide(); else show(t, true);
  });
  // keyboard focus (focus-visible only: a mouse press does not open it)
  root.addEventListener("focusin", e => {
    if (tip.contains(e.target)) { clearTimeout(hideTimer); return; }
    const t = trigOf(e);
    if (!t) return;
    if (quiet) return;
    let fv = false; try { fv = e.target.matches(":focus-visible"); } catch (_) { fv = true; }
    if (fv) show(t, false);
  });
  root.addEventListener("focusout", e => {
    if (!openTrig) return;
    const to = e.relatedTarget;
    if (to && (tip.contains(to) || openTrig.contains(to))) return;
    if (pinned && !(e.target === openTrig || openTrig.contains(e.target) || tip.contains(e.target))) return;
    hide();
  });
  // Esc closes; captured at window level so reveal.js does not also open its overview
  const onEsc = e => {
    if (e.key !== "Escape" || tip.hidden) return;
    const back = tip.contains(document.activeElement) ? openTrig : null;
    hide();
    // return focus to the trigger without re-opening its bubble (Esc means "closed")
    if (back && back.focus) { quiet = true; back.focus({ preventScroll: true }); quiet = false; }
    e.preventDefault(); e.stopPropagation(); e.stopImmediatePropagation();
  };
  window.addEventListener("keydown", onEsc, true);
  // a tap or click anywhere outside the widget closes a pinned bubble
  const onDocDown = e => { if (!tip.hidden && !root.contains(e.target)) hide(); };
  document.addEventListener("pointerdown", onDocDown, true);
  // stop listening once the widget leaves the page (OJS re-runs replace the node)
  const mo = new MutationObserver(() => { if (!root.isConnected && root._cx2mounted) { window.removeEventListener("keydown", onEsc, true); document.removeEventListener("pointerdown", onDocDown, true); mo.disconnect(); } if (root.isConnected) root._cx2mounted = true; });
  mo.observe(document.documentElement, { childList: true, subtree: true });

  // ---- wire controls -----------------------------------------------------------
  wC.s.addEventListener("input", () => { w = clean(snap(wC.s.value, G.w)); render(); });
  rC.s.addEventListener("input", () => { rho = clean(snap(rC.s.value, G.rho)); render(); });
  vC.s.addEventListener("input", () => { ratio = clean(snap(vC.s.value, G.ratio)); render(); });
  // lands exactly on w* (off the 0.01 grid on purpose); clampW is a no-op for every w* the
  // controls can produce, and a clamp, if the QA hook forces one, is flagged in readout + bubble
  toStarB.onclick = () => { w = clampW(K.wStar(V1, ratio * V1, rho)); render(); };
  // exact threshold (off the slider grid on purpose; the readout shows three decimals).
  // aria-disabled rather than disabled, so the button stays focusable and its bubble can say why.
  toThrB.onclick = () => { if (thrAvail()) { rho = K.rhoThreshold(V1, ratio * V1); render(); } };
  presetBtns.forEach(([b, p]) => { b.onclick = () => { rho = p.set.rho; ratio = p.set.ratio; render(); }; });
  resetB.onclick = () => { w = DEF.w; rho = DEF.rho; ratio = DEF.ratio; render(); };
  // Keys that reveal.js treats as navigation stop at the widget, so a focused slider or
  // button moves itself and not the slide.
  const NAV = ["ArrowLeft", "ArrowRight", "ArrowUp", "ArrowDown", " ", "Home", "End", "PageUp", "PageDown"];
  root.addEventListener("keydown", e => { if (NAV.indexOf(e.key) >= 0) e.stopPropagation(); });
  root.addEventListener("keyup", e => { if (NAV.indexOf(e.key) >= 0) e.stopPropagation(); });

  render();

  // hooks for headless QA (not part of the API)
  root._cx2 = { get state() { return { w, rho, ratio }; }, get regime() { return root.dataset.regime; },
                get checkFailed() { return checkFailed; }, get problems() { return problems.slice(); },
                get tipKeys() { return Object.keys(TIPS); }, get refs() { return JSON.parse(JSON.stringify(REFS)); },
                get openTip() { return tip.hidden ? null : { key: openKey, pinned, text: tip.textContent }; },
                show: key => { const t = root.querySelector(`[data-tip="${key}"]`); if (t) show(t, true); return !!t; }, hide,
                set(s) { if (s.w != null) w = s.w; if (s.rho != null) rho = s.rho; if (s.ratio != null) ratio = s.ratio; render(); } };
  return root;
}

// Payload: notes §7.7 Stages 2-3, per-replication holdout RMSE / simple-average RMSE.
// Everything between the markers is written by make_cpuzzle_data.R; do not edit by hand.
cpuzzleData = (
// CPUZZLE DATA BEGIN
{
  "source": "notes/module_07_notes.qmd section 7.7 Stages 2-3 (chunks stage23, stage3), via make_cpuzzle_data.R",
  "seed": 6378, "R": 200, "holdout": 12, "M": 6,
  "Tw": [12,24,60,120],
  "methods": ["IMSE", "CLS", "GR"],
  "arima1_median": [0.897666,0.901201,0.920388,0.897341],
  "summary": {
    "IMSE": {"median": [0.985218,0.984453,0.986148,0.982666], "p90": [1.014157,1.012969,1.011640,1.008623], "share": [0.745,0.780,0.780,0.830]},
    "CLS": {"median": [0.950174,0.928616,0.935850,0.900524], "p90": [1.120660,1.055176,1.050491,1.025199], "share": [0.640,0.730,0.755,0.855]},
    "GR": {"median": [1.677511,1.115839,0.984997,0.926358], "p90": [4.070284,1.922519,1.256734,1.108591], "share": [0.160,0.370,0.560,0.690]}
  },
  "rel": {
    "12": {
      "IMSE": [0.963,0.994,1.000,0.940,0.962,0.981,0.859,0.960,0.977,0.987,0.905,0.967,1.053,0.992,0.999,0.985,0.987,0.975,1.005,0.982,0.921,1.022,1.005,0.957,1.011,0.978,0.999,0.935,0.923,0.980,1.002,0.919,0.973,1.000,1.000,0.997,0.991,0.993,0.995,1.006,1.031,1.001,0.952,0.991,0.995,1.017,0.982,0.970,0.923,0.882,0.991,0.961,0.904,1.004,0.992,1.008,0.970,0.983,0.967,0.993,0.984,0.996,0.984,0.993,0.925,0.994,0.972,0.985,1.023,0.883,0.956,0.994,0.959,0.979,0.967,0.936,0.961,1.016,0.992,0.993,0.948,0.909,0.994,1.033,0.932,1.005,0.996,0.978,1.001,0.990,0.990,1.045,1.045,0.975,1.028,0.928,0.987,0.989,0.989,0.976,1.052,1.022,0.884,1.072,0.887,1.008,1.001,1.017,0.950,1.007,0.994,0.991,0.984,1.029,0.960,0.907,1.013,0.951,0.990,0.980,0.932,0.961,1.000,0.965,0.973,0.989,0.990,0.970,1.014,1.031,0.986,0.983,0.946,1.043,1.013,1.005,0.996,0.997,0.966,0.988,0.975,0.968,1.002,0.990,0.972,0.924,0.969,0.941,1.001,1.009,1.023,0.952,0.996,0.916,0.900,1.012,0.901,0.960,0.954,0.957,0.988,0.992,0.877,0.980,0.951,0.955,0.985,0.980,0.915,1.002,0.990,1.047,0.927,0.991,1.014,0.991,0.982,1.006,0.985,0.808,0.957,1.001,0.991,0.866,0.974,0.935,0.970,0.931,1.009,0.985,0.974,1.029,0.963,0.967,0.950,0.989,1.001,0.977,1.002,0.958],
      "CLS": [0.789,1.022,0.993,0.748,0.861,1.005,0.663,0.932,0.928,1.068,0.766,0.940,1.120,0.846,0.998,0.814,0.917,0.838,0.980,0.916,0.874,1.076,1.080,0.807,0.946,0.956,1.240,0.831,0.866,0.834,0.945,0.841,0.771,1.003,1.061,1.018,1.033,1.185,0.974,1.041,1.096,1.055,0.859,0.993,0.976,1.086,0.957,0.962,0.629,0.777,0.995,0.870,0.602,0.932,1.032,0.945,0.864,0.865,0.933,0.868,0.978,1.006,0.959,1.139,0.782,1.008,0.850,0.942,1.163,0.599,1.022,0.846,0.822,0.937,0.851,0.870,0.923,1.086,1.057,1.000,0.920,0.726,1.011,1.249,0.793,1.201,1.104,0.773,0.984,0.979,0.938,1.313,1.420,0.982,1.394,0.906,0.886,0.754,1.023,0.842,1.240,0.967,0.801,1.227,0.909,1.061,0.978,1.127,0.841,1.032,1.100,0.932,1.030,1.019,1.027,0.802,0.912,0.763,0.993,1.057,1.176,0.763,0.985,1.042,0.833,0.940,0.980,1.004,1.051,1.054,0.953,1.063,0.779,1.475,1.182,1.041,1.249,1.001,0.895,1.003,0.870,0.869,1.070,0.971,0.803,0.835,0.917,0.691,0.965,1.013,0.960,0.858,1.060,0.827,0.753,1.085,0.586,0.980,0.875,1.007,0.935,0.927,0.777,0.811,0.784,0.858,0.926,1.037,0.780,0.958,0.929,1.303,0.776,1.024,1.232,1.245,0.942,1.014,0.902,0.859,0.754,1.014,0.900,0.725,0.864,0.861,0.873,0.681,1.166,1.014,0.900,1.119,0.834,0.766,1.028,1.100,0.977,0.948,0.945,0.965],
      "GR": [0.922,3.271,1.899,1.208,1.195,3.438,3.285,1.791,0.967,3.328,0.961,2.202,3.570,6.633,2.080,0.660,0.973,3.011,1.980,1.795,4.106,1.926,1.677,2.455,0.997,1.440,2.917,3.177,0.943,10.595,0.905,1.638,1.089,2.201,1.516,1.479,1.409,1.808,0.944,4.121,1.772,5.516,3.867,1.297,1.865,1.459,1.920,1.293,0.925,0.832,1.415,2.447,1.349,1.089,1.996,1.085,5.782,1.043,1.472,0.989,3.547,1.127,1.283,1.652,0.933,1.213,4.066,2.329,2.835,1.575,1.007,2.812,2.135,1.678,0.736,1.295,0.916,1.865,1.313,1.136,3.452,1.674,2.650,4.636,0.772,2.759,4.164,1.963,1.244,1.568,1.433,2.179,5.520,1.399,2.675,1.289,2.060,1.787,3.540,0.950,1.145,2.120,0.938,1.648,1.621,2.266,1.314,3.776,0.670,2.003,11.683,3.754,3.252,2.807,1.947,0.768,0.947,1.376,0.920,1.576,1.526,1.602,2.450,1.500,0.932,1.231,5.741,2.894,1.982,1.738,0.937,0.899,1.766,4.828,3.185,1.280,3.795,0.969,0.999,1.324,1.806,3.318,2.737,1.088,1.104,0.955,4.160,1.989,1.041,2.476,1.047,5.291,2.319,1.071,0.963,3.117,2.110,1.093,1.410,1.004,1.068,1.968,0.884,0.922,1.549,1.215,2.336,8.028,2.704,1.508,1.953,6.889,1.061,1.330,2.368,5.915,1.111,1.104,1.678,2.877,3.614,2.610,2.428,1.300,0.932,3.244,1.625,4.259,4.372,2.294,1.157,1.379,1.830,1.476,1.119,10.017,1.352,1.849,2.670,2.098]
    },
    "24": {
      "IMSE": [0.987,0.983,1.015,0.980,0.906,0.964,1.022,1.000,0.969,1.011,0.927,0.999,0.990,0.993,0.983,0.980,0.984,0.972,0.987,1.029,0.931,0.996,1.012,0.994,0.990,0.994,0.985,1.026,0.906,0.972,0.974,0.936,0.987,0.985,1.004,1.012,0.995,0.978,0.935,0.999,0.981,0.999,0.941,0.968,0.944,1.004,0.976,1.007,0.955,1.100,0.952,0.957,0.939,0.992,0.973,0.998,0.981,0.963,1.031,1.005,1.001,0.975,1.009,1.007,0.950,0.993,0.981,0.995,1.007,0.947,0.982,0.992,0.963,0.973,0.930,0.941,1.009,0.938,0.981,0.970,0.970,1.033,0.996,1.008,0.962,0.978,0.967,0.955,1.013,0.993,0.992,0.935,1.006,0.978,1.004,0.997,0.912,0.957,0.995,0.989,1.012,0.998,0.960,1.019,1.029,0.987,0.998,0.978,1.017,0.980,0.993,1.030,0.937,1.011,0.970,1.040,0.991,0.910,0.999,0.987,0.990,0.970,0.998,0.947,0.989,1.004,1.004,1.066,0.982,1.020,0.999,0.976,0.943,0.972,0.887,0.998,0.985,0.983,0.948,1.010,0.997,0.990,0.994,0.976,0.998,0.975,1.003,0.929,0.994,0.968,1.018,0.960,0.979,0.955,0.864,0.989,0.985,0.949,0.976,1.032,0.959,0.964,1.037,0.915,0.959,0.929,0.996,0.938,0.973,1.013,0.964,0.993,0.956,0.985,0.990,0.990,0.981,0.999,0.964,0.894,0.961,0.985,0.998,0.976,0.989,0.964,0.962,1.016,0.980,0.976,0.982,1.006,0.969,0.918,0.981,1.011,0.949,0.951,0.959,1.038],
      "CLS": [0.894,0.992,0.897,1.073,0.681,0.891,1.036,0.850,0.814,1.004,0.968,1.081,0.941,1.009,1.002,0.903,0.921,0.754,0.878,1.054,0.943,0.981,1.064,1.023,0.940,0.981,1.054,0.918,0.656,0.826,0.856,0.785,1.045,0.929,1.054,1.001,1.007,0.847,0.920,0.871,0.971,0.956,1.004,0.811,0.872,0.959,0.844,1.043,0.844,1.247,0.668,0.819,0.814,1.028,0.897,0.995,0.920,0.753,0.978,0.933,1.000,0.760,1.001,1.053,0.934,0.881,0.889,0.957,1.020,0.921,0.953,1.002,0.821,0.858,0.692,0.780,1.012,0.898,1.038,0.920,0.835,1.081,1.141,1.090,0.931,0.873,0.944,0.708,1.035,0.976,1.026,0.806,1.012,0.861,0.999,1.046,0.913,0.712,1.014,1.064,0.986,0.907,0.877,0.898,1.004,0.847,0.978,0.985,0.898,0.928,0.941,1.210,0.798,0.926,0.873,1.082,0.916,0.639,0.974,0.951,0.983,0.871,1.051,0.860,0.880,0.920,0.976,1.099,0.914,1.143,1.104,0.872,0.848,1.034,0.765,1.120,1.053,0.830,0.755,0.995,1.192,1.030,0.847,0.877,0.950,0.927,0.846,0.778,1.047,0.992,1.067,0.956,0.860,0.882,0.662,0.976,0.885,0.794,0.810,1.090,0.836,0.721,1.040,0.654,0.913,0.948,0.964,0.826,0.968,1.010,0.908,0.989,0.832,0.896,0.884,0.940,0.870,0.963,0.937,1.008,0.832,0.971,1.265,0.898,0.979,0.928,0.905,1.087,0.904,0.947,0.857,0.914,0.851,0.891,0.973,1.092,0.848,0.838,0.714,0.927],
      "GR": [0.788,0.936,1.130,2.203,0.923,1.169,1.141,2.293,1.134,0.928,2.734,1.185,1.230,1.008,1.024,0.823,0.967,1.166,0.798,1.121,1.265,1.699,1.793,1.203,1.612,1.006,1.565,0.873,0.684,0.798,1.479,0.707,1.622,0.841,1.058,0.925,1.362,0.739,0.902,0.707,2.681,0.987,1.872,0.754,0.994,1.116,1.469,1.162,1.034,1.880,0.693,0.698,0.812,1.915,1.044,2.113,0.968,0.972,1.634,1.025,1.153,0.789,1.134,1.744,0.932,2.041,0.830,0.981,1.147,0.876,0.969,1.993,1.291,0.859,0.659,1.567,1.576,1.909,1.125,1.252,0.812,1.225,1.079,1.941,1.326,1.493,1.273,0.496,0.887,1.404,1.495,0.901,1.112,0.977,1.001,1.126,1.971,1.054,1.097,0.911,0.723,1.108,0.966,0.908,1.680,0.926,0.893,2.235,1.006,1.429,1.019,1.948,1.642,0.885,0.898,1.211,1.242,0.748,1.192,1.383,1.848,1.435,2.291,2.063,0.974,2.995,1.384,1.589,1.056,0.875,1.358,0.710,1.120,1.543,1.308,1.642,1.398,0.797,0.834,1.462,1.393,1.093,1.421,0.881,3.237,0.846,1.038,0.644,2.448,1.440,1.079,1.255,0.987,1.677,1.116,0.919,0.832,0.741,0.870,1.229,0.815,0.949,1.115,1.250,0.948,1.174,0.958,0.805,1.182,1.003,1.069,1.007,1.519,1.243,0.958,0.957,0.844,1.009,1.285,2.218,0.809,1.137,2.258,0.917,1.473,0.966,2.896,1.136,0.810,1.233,1.159,0.784,1.416,1.359,1.802,1.920,0.941,1.115,0.900,2.447]
    },
    "60": {
      "IMSE": [0.972,0.965,0.971,0.973,0.962,0.976,0.996,0.960,0.926,0.983,0.961,0.986,0.993,1.024,0.978,1.020,0.989,0.974,0.969,0.982,0.982,0.999,1.003,1.073,1.001,0.967,0.965,0.912,0.953,0.972,0.977,0.991,0.952,0.969,0.982,0.982,0.987,0.991,1.036,0.995,1.003,0.992,1.005,0.939,0.950,0.979,0.974,0.969,0.988,1.001,0.955,0.985,1.000,0.989,0.954,0.984,0.977,0.962,1.010,0.994,0.955,0.999,0.988,0.999,0.989,0.992,0.989,1.007,0.987,0.993,0.999,0.991,0.983,0.991,0.982,0.964,1.007,1.003,0.985,1.012,0.923,0.952,0.986,0.998,0.980,0.990,0.972,0.968,0.999,1.003,0.985,1.040,1.007,0.983,0.983,0.996,0.994,0.984,0.971,0.984,0.958,1.002,0.997,0.983,1.012,0.968,0.999,0.997,0.990,0.981,0.960,0.995,1.012,0.968,1.004,0.957,0.974,1.018,0.988,0.999,0.904,0.963,1.025,1.025,1.021,0.985,0.981,1.003,1.019,0.983,1.004,1.006,0.982,1.027,0.951,0.975,1.038,0.986,1.001,0.986,0.995,0.972,0.993,0.966,0.956,0.994,1.003,0.911,1.062,0.977,0.997,0.994,0.965,0.936,0.941,1.004,0.972,0.999,0.985,0.990,0.966,0.929,0.987,0.955,1.021,0.977,0.970,1.013,1.003,0.980,0.963,0.968,0.959,0.922,0.972,0.994,1.016,0.990,1.006,1.017,1.008,0.991,0.997,0.981,0.991,0.989,1.007,0.968,0.986,0.929,0.970,0.997,0.941,0.999,0.988,0.994,0.981,1.013,0.996,0.964],
      "CLS": [0.835,0.856,0.944,0.776,0.924,0.951,0.876,0.936,0.858,0.901,0.984,0.956,0.981,1.026,0.848,1.039,0.922,0.913,0.850,0.933,0.939,1.030,0.982,1.084,1.055,0.825,0.874,0.916,0.940,0.980,0.874,1.029,0.761,0.769,0.890,1.025,0.957,0.949,1.011,0.935,0.987,0.852,0.987,0.869,0.945,1.016,0.939,0.866,0.912,0.928,0.912,0.857,1.009,0.964,0.840,0.814,0.948,0.867,0.908,0.908,0.857,1.056,0.886,0.975,0.840,1.023,0.920,0.919,0.953,1.138,1.171,1.052,0.903,0.988,1.089,0.883,0.999,0.912,0.908,0.945,0.733,0.859,0.971,0.902,0.901,1.000,0.920,0.777,0.919,1.155,1.021,1.036,1.024,0.933,0.889,1.122,0.987,0.869,0.932,0.901,0.744,0.997,0.848,0.833,0.990,0.720,1.012,1.003,0.852,0.882,0.843,1.016,0.957,0.664,1.004,0.711,0.771,0.972,1.106,1.039,0.832,0.798,1.036,1.049,1.166,1.005,0.846,0.912,0.969,0.858,0.922,0.940,0.947,1.000,0.990,0.745,0.927,0.971,1.002,0.892,0.973,0.795,0.958,0.776,0.926,1.090,1.050,0.794,1.116,0.926,1.021,0.978,0.824,0.999,0.897,0.936,0.938,1.002,1.055,1.044,0.835,0.805,0.812,0.823,1.265,0.985,0.861,1.050,1.062,0.907,0.829,0.872,0.772,0.775,0.673,0.965,1.277,0.893,0.945,0.971,1.120,0.819,0.958,0.990,0.956,0.989,1.265,1.001,0.819,0.972,0.885,0.925,0.797,1.011,0.841,1.069,0.875,0.929,0.957,0.942],
      "GR": [0.707,0.925,1.428,0.772,0.994,0.925,0.870,1.166,0.910,0.855,0.905,0.853,1.010,1.385,0.838,1.482,0.959,1.108,0.986,0.943,0.986,1.090,0.890,2.442,1.422,0.835,0.872,0.967,1.055,1.249,0.869,1.120,1.010,1.212,1.198,1.024,0.963,1.014,1.054,0.917,1.083,0.824,0.908,0.921,1.299,0.993,1.072,0.924,0.911,0.883,0.981,0.815,1.352,1.016,0.818,0.772,1.055,1.084,0.919,1.052,1.031,1.100,0.851,1.035,0.943,0.952,0.945,0.860,0.977,1.394,1.172,1.014,1.057,0.798,1.428,0.962,0.989,0.940,0.886,1.101,0.749,0.801,1.036,0.955,0.808,1.117,0.949,0.775,0.929,1.193,1.032,1.277,0.988,0.963,0.974,1.347,0.957,0.845,1.262,0.984,1.005,0.969,0.866,0.839,1.117,1.044,1.070,1.068,1.111,0.904,0.787,0.974,0.813,0.589,0.989,0.579,0.956,1.038,1.884,0.981,0.913,0.905,1.113,1.194,1.257,1.200,0.878,0.843,1.051,0.996,0.978,0.889,1.014,1.113,1.019,0.736,0.936,1.152,0.984,1.284,0.980,0.752,0.989,0.729,0.895,1.155,1.149,0.669,1.238,0.939,1.050,1.084,0.704,1.013,0.875,1.042,1.013,1.065,1.348,1.175,0.781,1.123,0.829,0.820,1.297,0.975,0.852,1.036,1.110,1.258,0.753,0.903,0.713,0.800,1.041,0.951,1.460,0.996,0.874,1.047,2.000,0.782,1.145,0.926,1.058,1.017,1.373,0.928,0.812,0.989,0.934,0.934,1.069,1.095,0.800,1.188,0.895,0.923,1.045,0.992]
    },
    "120": {
      "IMSE": [0.970,0.953,0.994,0.966,0.980,0.993,0.950,0.981,0.979,0.960,0.989,0.973,1.031,0.967,0.949,0.979,0.922,0.957,0.952,0.952,1.008,1.000,0.969,0.949,1.006,0.959,0.986,0.993,0.961,0.985,0.981,0.984,0.957,0.986,0.974,0.959,0.958,0.996,0.965,0.957,0.983,0.975,0.933,0.968,1.015,0.991,0.982,1.012,0.949,1.000,0.975,0.993,0.989,1.005,0.996,0.993,0.971,0.983,0.997,0.975,1.013,0.991,0.948,0.967,0.992,0.984,0.989,0.950,0.966,0.994,0.935,1.039,0.979,0.983,0.992,0.982,1.000,1.014,0.964,0.986,1.001,0.952,0.978,1.020,0.995,1.005,0.991,0.950,1.001,0.995,1.003,0.984,0.982,0.991,0.976,0.952,0.956,0.963,0.997,0.977,0.953,0.996,1.010,0.985,0.943,1.018,0.994,0.926,0.980,0.964,0.968,0.977,1.001,0.990,0.966,0.981,1.000,0.925,0.984,0.989,1.011,0.998,0.976,0.988,0.987,1.018,0.956,0.988,1.008,0.979,0.978,1.026,0.943,0.972,0.998,0.957,0.969,0.975,0.948,0.994,0.993,0.988,0.986,1.014,0.983,0.964,0.966,0.972,1.033,0.990,1.006,0.993,1.010,0.987,0.956,0.971,0.972,0.954,0.975,1.001,0.961,0.990,0.963,0.990,0.990,1.011,0.963,0.995,0.968,0.969,0.981,0.972,0.956,1.023,0.959,0.938,1.001,0.985,0.988,0.995,0.986,0.993,0.923,0.987,0.963,0.982,0.985,0.983,0.991,0.945,1.016,0.918,0.989,0.964,1.019,0.985,0.993,1.015,0.974,0.935],
      "CLS": [0.901,0.822,0.886,0.840,0.799,1.095,0.990,0.806,0.803,0.803,1.013,0.957,0.990,0.878,0.931,0.921,0.776,0.945,0.818,0.714,0.983,1.030,0.751,0.887,0.959,0.829,0.939,1.208,0.828,0.876,0.821,0.927,0.694,0.889,0.983,0.977,0.814,1.033,0.840,0.717,0.896,0.948,0.817,0.890,0.938,0.875,0.801,0.941,1.012,0.986,0.858,0.965,0.953,1.061,0.903,0.991,0.761,0.908,0.930,0.943,0.947,0.894,0.695,0.897,1.010,0.833,0.939,0.982,0.778,0.926,0.776,1.094,0.828,0.951,1.003,0.900,0.924,0.970,1.045,0.946,1.025,0.810,0.934,1.102,1.007,0.900,0.921,0.747,0.920,0.908,1.136,0.897,0.831,0.924,0.913,0.839,0.818,0.809,0.775,0.878,0.682,0.896,0.979,0.892,0.924,1.048,0.963,0.946,0.884,0.784,0.788,0.839,0.905,0.997,0.876,0.884,0.954,0.852,0.857,0.963,0.974,0.874,0.869,0.877,0.904,0.941,0.809,0.896,1.158,0.773,0.873,1.062,0.858,0.923,0.776,0.815,0.941,0.991,0.845,0.975,0.824,0.922,0.810,1.100,0.878,0.798,0.696,0.924,0.978,0.958,0.944,0.984,1.105,0.953,0.875,1.004,0.771,0.725,0.838,1.049,0.592,0.879,0.849,0.986,0.833,0.965,0.797,0.945,0.870,0.776,0.910,0.796,0.832,1.006,0.741,0.931,1.037,0.970,0.916,0.838,0.789,0.863,0.907,0.807,0.782,0.895,0.837,0.885,1.003,0.858,1.095,0.709,0.898,0.993,1.105,0.936,0.991,1.048,1.076,0.678],
      "GR": [0.933,0.858,0.929,0.846,0.787,1.095,1.162,0.862,0.751,0.731,0.983,0.945,0.925,0.927,0.805,0.963,0.623,1.020,0.897,0.677,1.026,1.122,0.750,1.012,0.939,0.761,0.976,1.298,0.780,0.898,0.846,0.910,0.661,0.923,1.107,1.046,0.798,1.049,1.020,0.680,0.925,0.997,0.948,0.934,0.876,0.819,0.927,0.935,1.063,1.008,0.792,1.068,1.013,1.020,0.934,0.969,0.770,0.906,0.902,1.087,0.951,0.918,0.727,0.871,1.007,0.913,1.011,1.007,0.772,0.805,0.879,1.161,0.851,0.983,1.036,0.968,0.895,0.935,1.050,0.941,1.071,0.782,0.996,1.171,0.992,0.801,0.927,0.732,0.910,1.153,1.156,0.849,0.837,1.000,1.001,0.795,0.810,0.867,0.767,0.846,0.574,0.910,1.017,0.873,0.969,1.135,1.010,1.016,0.843,0.773,0.793,0.813,0.974,1.160,0.857,0.969,0.888,0.987,1.022,0.866,0.992,1.145,0.805,1.020,0.898,0.929,0.805,1.017,1.170,0.851,0.863,1.345,0.799,1.184,0.678,0.830,1.026,1.047,0.887,0.990,0.780,0.909,0.783,1.157,0.861,0.795,0.771,0.912,0.973,1.011,0.872,0.972,1.125,0.955,0.967,1.063,0.880,0.682,0.823,1.033,0.649,0.913,1.082,1.033,0.830,1.000,0.787,0.967,0.926,0.760,0.965,0.790,1.063,1.059,0.684,0.943,0.922,0.899,1.049,0.813,0.753,1.260,1.270,0.668,0.793,0.881,0.821,0.850,1.020,0.822,1.097,0.825,0.881,1.061,1.148,0.900,1.120,1.020,1.578,0.599]
    }
  },
  "win": {
    "12": {"IMSE": "11011111111101111101100101111101110111100011101111111010111111111111011111111011111010110110010111110010100010111011011111011111001110001111110111110001111011111111111110101101101110111111011011110101", "CLS": "10111011101101111111100111011111100000100011101111111101111110101011010111111000110010011110010111010110101010010001111001101110001010000010110111111011011011101111111011101000101110111111001011001111", "GR": "10000000101000011000000010001010000000100000000011000000000100001000000000101000000010000000000000010010000010000001101000001000001100000110000001000000001000000011000000000000000000001000000000000000"},
    "24": {"IMSE": "11011101101111111110110111101111110011111111101010111111110001001111011111110111101011110111010111110110011101101010111111111000101111111110111111011101111111101101111110111111111111111110111011101110", "CLS": "11101101101010011110110011011111010001111101111010111011111111001111011011110101100011110101011011001111011111101110111111011110100110100111001111110101111111101101111110111111111011011110111111101111", "GR": "11001000010000011010000000011101010101110101100000111000110001001011011001100000100000011001010000011011011000000110010000001000010100000110000101010000100111101100101100000011100010010100100100001010"},
    "60": {"IMSE": "11111111111110101111110001111111111111010101111110111111110111111110111111110010111111111010011111111011011111110101101111000110010010110101111111010111111011111111011001111111010001111101111111111011", "CLS": "11111111111110101111101001111110111011011111101111110111111110111011100011011111111111111000011011111111110011101101110011000011111110111101111110010101111110001111011001111111011101111100111110101111", "GR": "11011110111100101011101001110010000010010111010111110011001000101111100001011110110110111000111011010111000001111111100111000011011100011010111110010100101000001011011000111101011001010001111100101101"},
    "120": {"IMSE": "11111111111101111111001101111111111111111111011010111011111101111111111011111011011010110101111111111101101111110111111101111011011011111111111011110101011111101111101111111011011111111111110111011011", "CLS": "11111011110111111111101111101111111110111111111101111011111111110111111011011101011001111101111111111111101111111111111111111111011011111111111011111111011011101111111111111011011111111111010111011001", "GR": "11111001111111111011001011101111110010011111111100100011111011110100111011011101011011111001100111111101100011111011110110101110011010110011111011111011011011101100111111110011110110011111010110010001"}
  }
}
// CPUZZLE DATA END
)

// ============================== makeCombPuzzle ===============================
makeCombPuzzle = function (opts) {
  opts = opts || {};
  const D = opts.data || cpuzzleData;
  if (!D || !D.rel) throw new Error("makeCombPuzzle: payload missing; run make_cpuzzle_data.R");
  const uid = cpuzzleUid();

  // ---- parameters -----------------------------------------------------------
  const TW = D.Tw.slice();                                   // [12, 24, 60, 120]
  const MO = ["IMSE", "CLS", "GR"];                          // internal keys; never shown bare
  const NAME = { IMSE: "inverse-MSE", CLS: "constrained least squares", GR: "Granger–Ramanathan" };
  const LINES = { IMSE: ["inverse-MSE"], CLS: ["constrained", "least squares"], GR: ["Granger–", "Ramanathan"] };
  const KEY = { IMSE: "imse", CLS: "cls", GR: "gr" };
  const COL = { IMSE: "#2f8f6b", CLS: "#2c6fbb", GR: "#e4572e" };
  const DARK = { IMSE: "#1b5e45", CLS: "#1b4a82", GR: "#a8361a" };   // 90th-percentile bars
  const XMIN = opts.xmin == null ? 0.5 : +opts.xmin, XMAX = opts.xmax == null ? 5 : +opts.xmax;
  const reduce = typeof window !== "undefined" && window.matchMedia &&
                 window.matchMedia("(prefers-reduced-motion: reduce)").matches;
  const DUR = reduce ? 0 : (opts.duration == null ? 750 : Math.max(0, +opts.duration));
  let j = TW.indexOf(Number(opts.Tw0)); if (j < 0) j = 0;
  const j0 = j;
  let jShown = j;                                            // T_w the readouts describe
  const R = D.R;
  const val = (m, jj) => D.rel[String(TW[jj])][m];
  const S = (m, stat, jj) => D.summary[m][stat][jj];
  // full-precision "beats the simple average" flag (rel < 1 before the 3-dp rounding); a
  // payload without D.win (older make_cpuzzle_data.R) falls back to the rounded value
  const beat = (m, jj, i) => D.win ? D.win[String(TW[jj])][m].charAt(i) === "1" : val(m, jj)[i] < 1;

  // ---- "Defined on:" pointers --------------------------------------------------
  const STAGE2 = "Notes §7.7 — Stage 2, Monte Carlo (R = 200)";
  const REF0 = {
    Tw:         { label: "Notes §7.1 — The PS3 Design (the weights window)" },
    simple_avg: { label: "Notes §7.6 — Method 1, Simple Average" },
    imse:       { label: "Notes §7.6 — Method 2, Inverse-MSE Weights" },
    gr:         { label: "Notes §7.6 — Method 3, Granger-Ramanathan" },
    cls:        { label: "Notes §7.6 — Method 4, Constrained Least Squares" },
    rel_rmse:   { label: STAGE2 },
    median:     { label: STAGE2 },
    p90:        { label: STAGE2 },
    share:      { label: STAGE2 },
    clip:       { label: "Notes §7.7 — Stage 2 (Granger-Ramanathan's long right tail)" },
    true_model: { label: "Notes §7.1 — The Course Simulation (arima_1, the correctly specified model)" }
  };
  const REFS = {};
  Object.keys(REF0).forEach(k => { REFS[k] = Object.assign({}, REF0[k], (opts.refs && opts.refs[k]) || {}); });

  // ---- deterministic jitter (the course generator, verbatim from the M5/M6 cores) ----
  const mulberry32 = a => () => {
    a |= 0; a = a + 0x6D2B79F5 | 0;
    let t = Math.imul(a ^ a >>> 15, 1 | a);
    t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t;
    return ((t ^ t >>> 14) >>> 0) / 4294967296;
  };
  const hashSeed = a => {
    a = a >>> 0;
    a = Math.imul(a ^ (a >>> 16), 2246822507);
    a = Math.imul(a ^ (a >>> 13), 3266489909);
    return (a ^ (a >>> 16)) >>> 0;
  };
  const rng = mulberry32(hashSeed(D.seed));
  const JY = new Array(R), JX = new Array(R);                // one pair per replication
  for (let r = 0; r < R; r++) { JY[r] = 2 * rng() - 1; JX[r] = 2 * rng() - 1; }

  // ---- formatting and notation ------------------------------------------------
  const f2 = x => Number(x).toFixed(2);
  // Shares are k / 200, so 100 x share has at most one decimal (.5): print it exactly
  // ("74.5%", "16%") rather than rounding; R's round-half-even prints 0.745 as 0.74 in the
  // notes' Stage 2 table, JS's Math.round would print 75: exact sidesteps the conflict.
  const pct = x => { const v = Math.round(1000 * x) / 10; return (Number.isInteger(v) ? String(v) : v.toFixed(1)) + "%"; };
  const esc = s => String(s).replace(/&/g, "&amp;").replace(/</g, "&lt;").replace(/>/g, "&gt;").replace(/"/g, "&quot;");
  const M_ = s => `<math>${s}</math>`;
  const TwM = M_(`<msub><mi>T</mi><mi>w</mi></msub>`);
  const hatw = sub => M_(`<msub><mover accent="true"><mi>w</mi><mo>^</mo></mover>${sub}</msub>`);
  const wm = M_(`<msub><mi>w</mi><mi>m</mi></msub>`);
  const yt = M_(`<msub><mi>y</mi><mi>t</mi></msub>`);
  const Mm = M_(`<mi>M</mi>`);

  // ---- bubble texts: { title, body, value? } in course notation ----------------------
  function tipContent(key, el) {
    const tw = TW[jShown], m = el && el.dataset ? el.dataset.m : null;
    const rowVal = (stat, f) => m ? `<span class="cpuzzle-tv">${esc(NAME[m])} at ${TwM} = ${tw}: <b>${f(S(m, stat, jShown))}</b></span>` : "";
    switch (key) {
      case "Tw": {
        const v = el && el.dataset.tw ? +el.dataset.tw : tw;
        return { ref: "Tw", title: `Weights window, ${TwM} = ${v}`,
          body: `The ${v} dates on which the combination weights ${hatw(M_("<mi>m</mi>"))} are estimated; every method is then scored on the next 12 holdout dates. Problem Set 3 uses ${TwM} = 12.` };
      }
      case "simple_avg":
        return { ref: "simple_avg", title: "Simple average (SA)",
          body: `Every forecaster gets ${wm} = 1/${Mm}, here ${Mm} = ${D.M} tracks, and no intercept. Nothing is estimated, so it carries no estimation error: the benchmark at 1 on this axis.` };
      case "imse":
        return { ref: "imse", title: "Inverse-MSE weights (IMSE)",
          body: `${hatw(M_("<mi>m</mi>"))} ∝ 1/MSE<sub>m</sub>, each track's own weights-window mean squared error; non-negative, summing to one, no intercept. It estimates ${Mm} numbers, one column at a time.` };
      case "cls":
        return { ref: "cls", title: "Constrained least squares (CLS)",
          body: `Least squares of ${yt} on the ${Mm} forecasts with no intercept, ${hatw(M_("<mi>m</mi>"))} ≥ 0 and Σ ${hatw(M_("<mi>m</mi>"))} = 1. It estimates ${Mm} weights jointly, boxed in.` };
      case "gr":
        return { ref: "gr", title: "Granger–Ramanathan (GR)",
          body: `Regress ${yt} on all ${Mm} forecasts with an intercept and no constraints: ${Mm} + 1 numbers (${hatw(M_("<mn>0</mn>"))}, ${hatw(M_("<mn>1</mn>"))}, …, ${hatw(M_("<mi>M</mi>"))}) from ${TwM} dates. It fits the weights window best and carries the most estimation error.` };
      case "rel_rmse": {
        let v = "";
        if (el && el.dataset.rep && m) {
          // Classify with the full-precision flag, never the 3-dp value: 16 replications print
          // 1.000, half of them wins (e.g. 0.999848) and half losses.
          const r = +el.dataset.rep, x = val(m, jShown)[r - 1], w3 = x.toFixed(3);
          const won = beat(m, jShown, r - 1);
          const verdict = w3 === "1.000"
            ? (won ? "(beats the simple average at full precision)" : "(loses to the simple average at full precision)")
            : (won ? "(beat the simple average)" : "(lost to the simple average)");
          v = `<span class="cpuzzle-tv" data-win="${won ? 1 : 0}">Replication ${r}, ${esc(NAME[m])}, ${TwM} = ${tw}: <b>${w3}×</b> ${verdict}</span>`;
        }
        return { ref: "rel_rmse", title: "One replication",
          body: `One fresh series from the course DGP. The dot sits at the method's holdout RMSE divided by the simple average's holdout RMSE on the same 12 dates; below 1, the method beat the simple average. The same ${R} series at every ${TwM}.`, value: v };
      }
      case "beats":
        return { ref: "rel_rmse", title: "Beats the simple average",
          body: `Left of the dashed line the ratio is below 1: in that replication the combination's holdout RMSE was smaller than the simple average's.` };
      case "median":
        return { ref: "median", title: "Median (navy bar)",
          body: `The middle of the ${R} replications' ratios: half the series did better than this, half worse.`, value: rowVal("median", x => f2(x) + "×") };
      case "p90":
        return { ref: "p90", title: "90th percentile (dark coloured bar)",
          body: `One replication in ten did worse than this. Granger–Ramanathan's long right tail at short ${TwM} is the combination puzzle in one picture.`, value: rowVal("p90", x => f2(x) + "×") };
      case "share":
        return { ref: "share", title: "Share beating the simple average",
          body: `The fraction of the ${R} replications with a ratio below 1: the method's holdout RMSE was smaller than the simple average's.`, value: rowVal("share", pct) };
      case "clip": {
        let v = "";
        if (m) {
          const xs = val(m, jShown), over = xs.filter(x => x > XMAX), under = xs.filter(x => x < XMIN);
          v = `<span class="cpuzzle-tv">${esc(NAME[m])} at ${TwM} = ${tw}: <b>${over.length}</b> above ${XMAX}×${over.length ? ` (largest ${Math.max(...over).toFixed(2)}×)` : ""}, <b>${under.length}</b> below ${XMIN}×.</span>`;
        }
        return { ref: "clip", title: "Beyond the axis",
          body: `Ratios above ${XMAX}× are piled in the hatched lane with their count (below ${XMIN}×, pinned at the left edge). They are clipped only in the picture: the median, 90th percentile and share use every replication.`, value: v };
      }
      case "true_model":
        return { ref: "true_model", title: "Correctly specified model, arima_1 (grey dotted line)",
          body: `The median ratio for arima_1 used alone, with no combining: the correctly specified ARIMA(1,0,0)(1,0,0)<sub>12</sub> form, its coefficients estimated on each series like every other track. One reference value, drawn across all three rows.`,
          value: `<span class="cpuzzle-tv">At ${TwM} = ${tw}: <b>${f2(D.arima1_median[jShown])}×</b> the simple average</span>` };
    }
    return null;
  }

  // ---- DOM ------------------------------------------------------------------
  const el = (tag, cls, html) => { const n = document.createElement(tag); if (cls) n.className = cls; if (html != null) n.innerHTML = html; return n; };
  const root = el("div", "cpuzzle");
  root.dataset.widget = "cpuzzle"; root.dataset.uid = String(uid);

  const bar = el("div", "cpuzzle-bar");
  const grp = el("div", "cpuzzle-grp");
  grp.setAttribute("role", "radiogroup"); grp.setAttribute("aria-label", "weights-window length T_w");
  const BTN = TW.map((tw, k) => {
    const b = document.createElement("button");
    b.type = "button"; b.className = "cpuzzle-btn"; b.textContent = String(tw);
    b.dataset.action = "tw-" + tw; b.dataset.tip = "Tw"; b.dataset.tw = String(tw);
    b.setAttribute("role", "radio"); b.id = `cpuzzle${uid}-tw${tw}`;
    b.addEventListener("click", () => setJ(k));
    grp.appendChild(b); return b;
  });
  const lab = el("span", "cpuzzle-lab", `weights window ${TwM} =`);
  const tag = el("span", "cpuzzle-tag", `R = ${R} replications · holdout N = ${D.holdout} · seed ${D.seed}`);
  bar.append(lab, grp, tag);

  const chart = el("div", "cpuzzle-chart");
  // Legend: each item is a focusable bubble trigger (the SVG marks it names also open it on hover).
  const legend = el("div", "cpuzzle-legend");
  const sw = inner => `<svg viewBox="0 0 32 18" aria-hidden="true">${inner}</svg>`;
  const II = `<span class="cpuzzle-ii" aria-hidden="true">i</span>`;
  const LEG = [
    ["median", sw(`<line x1="16" y1="1" x2="16" y2="17" stroke="#0b1f3a" stroke-width="4" stroke-linecap="round"/>`), "median"],
    ["p90", sw(`<line x1="16" y1="1" x2="16" y2="17" stroke="#a8361a" stroke-width="3" stroke-linecap="round"/>`), "90th percentile"],
    ["rel_rmse", sw(`<circle cx="16" cy="9" r="5" fill="#2c6fbb" fill-opacity=".42" stroke="#2c6fbb"/>`), "one replication"],
    ["true_model", sw(`<line x1="16" y1="1" x2="16" y2="17" stroke="#6b7785" stroke-width="2.4" stroke-dasharray="2 3"/>`), "correctly specified (arima_1)"],
    ["beats", sw(`<rect x="2" y="1" width="28" height="16" fill="#2f8f6b" fill-opacity=".18"/>`), "beats the simple average"],
    ["clip", sw(`<rect x="2" y="1" width="28" height="16" fill="#fff" stroke="#b9c3cd"/><line x1="6" y1="17" x2="16" y2="1" stroke="#b9c3cd" stroke-width="2"/><line x1="16" y1="17" x2="26" y2="1" stroke="#b9c3cd" stroke-width="2"/>`), `beyond ${XMAX}×`]
  ];
  LEG.forEach(([k, s, t]) => {
    const it = el("span", "cpuzzle-li", `${s}<span>${t}</span>${II}`);
    it.tabIndex = 0; it.setAttribute("role", "button"); it.dataset.tip = k;
    it.setAttribute("aria-label", `About: ${t}`);
    legend.appendChild(it);
  });
  const cap = el("div", "cpuzzle-cap"); cap.dataset.role = "caption";
  const foot = el("div", "cpuzzle-foot");
  foot.innerHTML =
    `Each dot is one simulated series from the course DGP, the same ${R} series at every ${TwM}, so the dots slide as ` +
    `the window grows; the axis is on a log scale. Hover, focus or tap any item marked ${II} for a definition. Notes §7.7, Stages 2–3.`;
  const tip = el("div", "cpuzzle-tip"); tip.hidden = true; tip.id = `cpuzzle${uid}-tip`;
  tip.setAttribute("role", "tooltip");
  root.append(bar, chart, legend);
  if (opts.caption !== false) root.append(cap);
  if (opts.foot !== false) root.append(foot);
  root.append(tip);

  // ---- layout (rebuilt when the mount width or host font size changes) ----------
  const NS = "http://www.w3.org/2000/svg";
  let G = null;            // geometry
  let dots = null;         // dots[m][r] -> <circle>
  let parts = null;        // per-row markers and text
  let cur = MO.map(m => val(m, j).map(Math.log));     // displayed log values
  const lnMin = Math.log(XMIN), lnMax = Math.log(XMAX);
  const inDeck = () => !!root.closest(".reveal");

  function geometry() {
    const w = root.clientWidth;
    const W = Math.max(520, Math.round(w > 0 ? w : 1800));
    const host = parseFloat(getComputedStyle(root).fontSize) || 16;
    const deck = inDeck();
    // deck: 26px host -> fs 28.6; smallest SVG text 0.86 fs -> 24.6px CSS, ~24px projected.
    // notes: fs = the body size, so the narrow column keeps a usable plot width.
    const fs = host * (deck ? 1.1 : 1.0);
    const labelW = 8.6 * fs, statsW = 10.6 * fs, lane = 2.8 * fs, gap = 0.6 * fs;
    const xL = labelW, xR = W - statsW - lane - gap;
    const laneL = xR + 0.5 * fs, laneR = xR + lane;
    const top = 3.0 * fs, rowH = 4.6 * fs, bot = top + MO.length * rowH;
    const H = Math.round(bot + 3.0 * fs);
    const sx = lv => xL + (lv - lnMin) / (lnMax - lnMin) * (xR - xL);
    const statsX0 = W - statsW + 0.2 * fs;
    const colX = [statsX0 + 1.7 * fs, statsX0 + 5.2 * fs, statsX0 + 8.7 * fs];
    return { W, H, fs, xL, xR, laneL, laneR, top, rowH, bot, statsX0, colX, sx, host, w, deck };
  }

  function build() {
    G = geometry();
    const { W, H, fs, xL, xR, laneL, laneR, top, rowH, bot, sx, colX } = G;
    const F = v => v.toFixed(1);
    const T = (x, y, s, a, extra) => `<text x="${F(x)}" y="${F(y)}" text-anchor="${a || "start"}" ${extra || ""}>${s}</text>`;
    // focusable info mark: a small circled i
    const I = (cx, cy, key, label, extra) =>
      `<g class="cpuzzle-i" tabindex="0" role="button" aria-label="About: ${esc(label)}" data-tip="${key}" ${extra || ""} transform="translate(${F(cx)},${F(cy)})">` +
      `<circle r="${F(0.5 * fs)}" fill="#fff" stroke="#5d6b78" stroke-width="1.6"/>` +
      `<text y="${F(0.28 * fs)}" text-anchor="middle" font-size="${F(0.78 * fs)}" font-style="italic" font-weight="700" font-family="Georgia, serif" fill="#5d6b78" class="cpuzzle-glyph">i</text></g>`;
    let s = "";
    s += `<defs><pattern id="cpuzzle${uid}-hatch" width="${F(0.5 * fs)}" height="${F(0.5 * fs)}" patternUnits="userSpaceOnUse" patternTransform="rotate(45)">` +
         `<line x1="0" y1="0" x2="0" y2="${F(0.5 * fs)}" stroke="#d7dee5" stroke-width="${F(0.12 * fs)}"/></pattern></defs>`;
    // row bands, then the beats-SA tint (a bubble target)
    MO.forEach((m, i) => {
      const y0 = top + i * rowH;
      if (i % 2 === 1) s += `<rect x="0" y="${F(y0)}" width="${W}" height="${F(rowH)}" fill="#eef2f6" fill-opacity=".55"/>`;
    });
    const x1 = sx(0);
    s += `<rect data-tip="beats" x="${xL}" y="${top}" width="${F(x1 - xL)}" height="${F(bot - top)}" fill="#2f8f6b" fill-opacity=".07"/>`;
    // row labels: name (spelled out, one or two lines) with an info mark to its left
    MO.forEach((m, i) => {
      const y0 = top + i * rowH, L = LINES[m], lh = 1.12 * fs;
      const yFirst = y0 + rowH / 2 - (L.length - 1) * lh / 2 + 0.33 * fs;
      s += `<g data-tip="${KEY[m]}" data-k="label" data-m="${m}">`;
      L.forEach((t, q) => { s += T(1.35 * fs, yFirst + q * lh, esc(t), "start", `font-size="${F(0.98 * fs)}" font-weight="700" fill="${COL[m]}"`); });
      s += `</g>`;
      s += I(0.6 * fs, y0 + rowH / 2, KEY[m], NAME[m]);
    });
    // overflow lane (a bubble target)
    s += `<rect data-tip="clip" x="${F(laneL)}" y="${top}" width="${F(laneR - laneL)}" height="${F(bot - top)}" fill="url(#cpuzzle${uid}-hatch)" stroke="#d7dee5"/>`;
    s += T((laneL + laneR) / 2, top - 0.55 * fs, `&gt; ${XMAX}×`, "middle", `data-tip="clip" font-size="${F(0.86 * fs)}" fill="#5d6b78"`);
    // x ticks; narrow plot (notes column): drop the 0.75 and 1.5 ticks so the labels do not collide
    const narrow = (xR - xL) < 30 * fs;
    const ticks = (narrow ? [0.5, 1, 2, 3, 5] : [0.5, 0.75, 1, 1.5, 2, 3, 5]).filter(v => v >= XMIN - 1e-9 && v <= XMAX + 1e-9);
    ticks.forEach(v => {
      const x = sx(Math.log(v));
      s += `<line x1="${F(x)}" y1="${top}" x2="${F(x)}" y2="${bot}" stroke="#e7ecf1" stroke-width="1" pointer-events="none"/>`;
      s += `<line x1="${F(x)}" y1="${bot}" x2="${F(x)}" y2="${F(bot + 0.35 * fs)}" stroke="#8a94a0"/>`;
      s += T(x, bot + 1.35 * fs, (v === 1 ? "1" : String(v)) + "×", "middle", `font-size="${F(fs)}" fill="${v === 1 ? "#0b1f3a" : "#5d6b78"}" font-weight="${v === 1 ? 700 : 400}"`);
    });
    s += `<line x1="${xL}" y1="${bot}" x2="${xR}" y2="${bot}" stroke="#16202b" stroke-width="1.4"/>`;
    // direction words and axis title
    s += T(x1 - 0.5 * fs, bot + 2.6 * fs, "◂ beats the simple average", "end", `data-tip="beats" font-size="${F(0.88 * fs)}" fill="#2f8f6b"`);
    s += T(x1 + 0.5 * fs, bot + 2.6 * fs, "loses to it ▸", "start", `font-size="${F(0.88 * fs)}" fill="#8a94a0"`);
    s += T(colX[0] - 1.5 * fs, bot + 1.35 * fs, "holdout RMSE ÷", "start", `data-tip="rel_rmse" font-size="${F(0.86 * fs)}" fill="#5d6b78"`);
    s += T(colX[0] - 1.5 * fs, bot + 2.45 * fs, "simple-average RMSE", "start", `data-tip="rel_rmse" font-size="${F(0.86 * fs)}" fill="#5d6b78"`);
    // stats column headers, each with an info mark above it
    [["median", "median"], ["p90", "90th pct"], ["share", "beats SA"]].forEach(([k, t], c) => {
      s += T(colX[c], top - 0.55 * fs, t, "middle", `data-tip="${k}" data-k="h-${k}" font-size="${F(0.86 * fs)}" fill="#5d6b78" font-weight="600"`);
      s += I(colX[c], top - 2.3 * fs, k, t === "beats SA" ? "share beating the simple average" : t === "90th pct" ? "90th percentile" : "median");
    });
    // correctly-specified-model (arima_1) reference: one faint dotted segment per row, moved as a group
    s += `<g data-tip="true_model" data-k="tm">`;
    MO.forEach((m, i) => {
      const y0 = top + i * rowH;
      s += `<line x1="0" x2="0" y1="${F(y0 + 0.05 * rowH)}" y2="${F(y0 + 0.95 * rowH)}" stroke="#6b7785" stroke-width="${F(0.1 * fs)}" stroke-dasharray="${F(0.12 * fs)} ${F(0.22 * fs)}" stroke-opacity=".8"/>`;
    });
    s += `<line x1="0" x2="0" y1="${top}" y2="${bot}" stroke="transparent" stroke-width="${F(0.5 * fs)}"/></g>`;
    // per-row dynamic layers: dots group, markers, counts, stats
    MO.forEach((m, i) => {
      const y0 = top + i * rowH;
      s += `<g data-row="${m}">`;
      s += `<g class="cpuzzle-dots" data-method="${m}"></g>`;
      const bar = (k, y1, y2, stroke, wid) =>
        `<g data-tip="${k}" data-k="${k}" data-m="${m}">` +
        `<line x1="0" x2="0" y1="${F(y1)}" y2="${F(y2)}" stroke="${stroke}" stroke-width="${F(wid)}" stroke-linecap="round"/>` +
        `<line class="cpuzzle-hit" x1="0" x2="0" y1="${F(y1)}" y2="${F(y2)}" stroke="transparent" stroke-width="${F(0.7 * fs)}"/></g>`;
      s += bar("p90", y0 + 0.14 * rowH, y0 + 0.86 * rowH, DARK[m], 0.17 * fs);
      s += bar("median", y0 + 0.08 * rowH, y0 + 0.92 * rowH, "#0b1f3a", 0.24 * fs);
      s += T((laneL + laneR) / 2, y0 + 0.95 * fs, "", "middle", `data-tip="clip" data-m="${m}" data-k="over" font-size="${F(0.86 * fs)}" font-weight="700" fill="#c2401b"`);
      s += T(xL + 0.15 * fs, y0 + 0.95 * fs, "", "start", `data-tip="clip" data-m="${m}" data-k="under" font-size="${F(0.86 * fs)}" fill="#5d6b78"`);
      const yv = y0 + rowH / 2 + 0.33 * fs;
      s += T(colX[0], yv, "", "middle", `data-tip="median" data-m="${m}" data-k="s-med" font-size="${F(0.95 * fs)}" fill="#0b1f3a" font-weight="700"`);
      s += T(colX[1], yv, "", "middle", `data-tip="p90" data-m="${m}" data-k="s-p90" font-size="${F(0.95 * fs)}" fill="#16202b"`);
      s += T(colX[2], yv, "", "middle", `data-tip="share" data-m="${m}" data-k="s-shr" font-size="${F(0.95 * fs)}" fill="#16202b"`);
      s += `</g>`;
    });
    // the simple-average line (a bubble target, with a wide invisible hit stroke) and its label + mark
    s += `<g data-tip="simple_avg" data-k="sa">` +
         `<line class="cpuzzle-hit" x1="${F(x1)}" y1="${F(top - 0.2 * fs)}" x2="${F(x1)}" y2="${bot}" stroke="#0b1f3a" stroke-width="${F(0.12 * fs)}" stroke-dasharray="${F(0.45 * fs)} ${F(0.3 * fs)}"/>` +
         `<line class="cpuzzle-hit" x1="${F(x1)}" y1="${F(top - 0.2 * fs)}" x2="${F(x1)}" y2="${bot}" stroke="transparent" stroke-width="${F(0.45 * fs)}"/>` +
         T(x1, top - 0.55 * fs, "simple average", "middle", `font-size="${F(0.9 * fs)}" font-weight="700" fill="#0b1f3a"`) + `</g>`;
    s += I(x1 + 4.3 * fs, top - 0.85 * fs, "simple_avg", "simple average", `data-k="i-sa"`);

    chart.innerHTML = `<svg class="cpuzzle-svg" viewBox="0 0 ${W} ${H}" role="img" aria-label="Holdout RMSE relative to the simple average for inverse-MSE, constrained least squares and Granger-Ramanathan combinations, ${R} replications, log scale"><title>Combination puzzle: holdout RMSE relative to the simple average</title>${s}</svg>`;
    const svg = chart.querySelector("svg");
    dots = {}; parts = {};
    MO.forEach((m, i) => {
      const g = svg.querySelector(`g.cpuzzle-dots[data-method="${m}"]`);
      const row = svg.querySelector(`g[data-row="${m}"]`);
      dots[m] = [];
      for (let r = 0; r < R; r++) {
        const c = document.createElementNS(NS, "circle");
        c.setAttribute("r", (0.27 * fs).toFixed(1));
        c.setAttribute("fill", COL[m]); c.setAttribute("fill-opacity", ".42");
        c.setAttribute("stroke", COL[m]); c.setAttribute("stroke-opacity", ".85"); c.setAttribute("stroke-width", "1");
        c.dataset.rep = String(r + 1); c.dataset.m = m; c.dataset.tip = "rel_rmse";
        g.appendChild(c); dots[m].push(c);
      }
      const q = k => row.querySelector(`[data-k="${k}"]`);
      parts[m] = { med: q("median"), p90: q("p90"), over: q("over"), under: q("under"),
                   sMed: q("s-med"), sP90: q("s-p90"), sShr: q("s-shr"), y0: top + i * rowH };
    });
    parts.tm = svg.querySelector('[data-k="tm"]');
    G.lastW = root.clientWidth; G.lastHost = G.host;
  }

  // ---- draw the current (possibly mid-tween) state ---------------------------------
  // cur: displayed log ratios; mk: displayed log median / p90 per method, plus tm (arima_1, correctly specified)
  function draw(mk) {
    const { fs, xL, laneL, laneR, rowH, sx } = G;
    const clampX = lv => sx(Math.min(lnMax, Math.max(lnMin, lv)));
    MO.forEach((m, i) => {
      const P = parts[m], y0 = P.y0, yc = y0 + rowH / 2, amp = 0.36 * rowH;
      let over = 0, under = 0;
      const lv = cur[i];
      for (let r = 0; r < R; r++) {
        let x, y = yc + 0.12 * fs + JY[r] * amp * 0.82;
        if (lv[r] > lnMax) {                       // overflow lane: below the row's count
          over++; x = (laneL + laneR) / 2 + JX[r] * 0.3 * (laneR - laneL);
          y = y0 + 0.5 * rowH + (JY[r] + 1) / 2 * 0.38 * rowH;
        }
        else if (lv[r] < lnMin) { under++; x = xL + 0.3 * fs; }
        else x = sx(lv[r]);
        const c = dots[m][r];
        c.setAttribute("cx", x.toFixed(1));
        c.setAttribute("cy", y.toFixed(1));
      }
      P.med.setAttribute("transform", `translate(${clampX(mk.rows[i].med).toFixed(1)},0)`);
      P.p90.setAttribute("transform", `translate(${clampX(mk.rows[i].p90).toFixed(1)},0)`);
      P.over.textContent = over ? String(over) : "";
      P.over.dataset.n = String(over);
      P.under.textContent = under ? `◂ ${under} below ${XMIN}×` : "";
      P.under.dataset.n = String(under);
    });
    parts.tm.setAttribute("transform", `translate(${clampX(mk.tm).toFixed(1)},0)`);
  }

  // Readouts and caption describe jShown; called when a transition ENDS (or immediately when
  // there is no transition), so the numbers never disagree with the dots on screen.
  function words() {
    const jj = jShown;
    MO.forEach(m => {
      const P = parts[m];
      P.sMed.textContent = `${f2(S(m, "median", jj))}×`;
      P.sP90.textContent = `${f2(S(m, "p90", jj))}×`;
      P.sShr.textContent = pct(S(m, "share", jj));
    });
    const tw = TW[jj];
    const g = S("GR", "median", jj), gp = S("GR", "p90", jj), gs = S("GR", "share", jj);
    const c = S("CLS", "median", jj), im = S("IMSE", "median", jj);
    const GRs = s => `<span class="cpuzzle-gr">${s}</span>`;
    let line;
    // Fixed wording per T_w, numbers from the payload; each sentence is guarded by the
    // inequality it asserts, with a neutral fallback if a regenerated payload breaks it.
    if (tw === 12 && g > 1 && c < 1 && im < 1)
      line = `${GRs("Granger–Ramanathan's estimated weights overfit")}: at the median its holdout RMSE is ${f2(g)}× the simple average's (90th percentile ${f2(gp)}×). At the median the restricted methods already beat the simple average: constrained least squares ${f2(c)}×, inverse-MSE ${f2(im)}×.`;
    else if (tw === 24 && g > 1)
      line = `Twice the window and ${GRs("Granger–Ramanathan still loses")}: median ${f2(g)}× the simple average, beating it in only ${pct(gs)} of replications.`;
    else if (tw === 60 && g < 1)
      line = `Granger–Ramanathan crosses the simple average: median ${f2(g)}×, but its 90th percentile is still ${f2(gp)}×.`;
    else if (tw === 120 && g < 1)
      line = `Ten years of window and estimating weights pays: Granger–Ramanathan ${f2(g)}×, constrained least squares ${f2(c)}× (correctly specified arima_1 alone: ${f2(D.arima1_median[jj])}×).`;
    else
      line = `Medians relative to the simple average: Granger–Ramanathan ${f2(g)}×, constrained least squares ${f2(c)}×, inverse-MSE ${f2(im)}×.`;
    cap.innerHTML = `${TwM} = ${tw}: ${line}`;
    root.dataset.shown = String(tw);
    if (tipState.open) refreshTip();
  }
  function buttons() {
    BTN.forEach((b, k) => {
      const on = k === j;
      b.classList.toggle("cpuzzle-on", on); b.setAttribute("aria-checked", on ? "true" : "false"); b.tabIndex = on ? 0 : -1;
    });
    root.dataset.tw = String(TW[j]);
  }

  const markers = jj => ({ rows: MO.map(m => ({ med: Math.log(S(m, "median", jj)), p90: Math.log(S(m, "p90", jj)) })),
                           tm: Math.log(D.arima1_median[jj]) });
  let mkCur = markers(j);

  // ---- transitions -------------------------------------------------------------
  let raf = null;
  function setJ(k) {
    if (k < 0 || k >= TW.length || (k === j && !raf)) return;
    const from = cur.map(a => a.slice()), mkFrom = mkCur;
    j = k; buttons();
    const to = MO.map(m => val(m, j).map(Math.log)), mkTo = markers(j);
    if (raf) { cancelAnimationFrame(raf); raf = null; }
    const finish = () => { cur = to; mkCur = mkTo; draw(mkCur); jShown = j; words(); };
    if (DUR === 0) { finish(); return; }
    const t0 = performance.now();
    const ease = u => u < 0.5 ? 4 * u * u * u : 1 - Math.pow(-2 * u + 2, 3) / 2;
    const lerp = (a, b, e) => a + (b - a) * e;
    const step = now => {
      if (!root.isConnected) { raf = null; return; }          // detached: stop
      const u = Math.min(1, (now - t0) / DUR), e = ease(u);
      if (u >= 1) { raf = null; finish(); return; }
      cur = from.map((a, i) => a.map((v, r) => lerp(v, to[i][r], e)));
      mkCur = { rows: mkFrom.rows.map((o, i) => ({ med: lerp(o.med, mkTo.rows[i].med, e), p90: lerp(o.p90, mkTo.rows[i].p90, e) })),
                tm: lerp(mkFrom.tm, mkTo.tm, e) };
      draw(mkCur);
      raf = requestAnimationFrame(step);
    };
    raf = requestAnimationFrame(step);
  }

  // ---- info bubbles --------------------------------------------------------------
  const tipState = { open: false, pinned: false, target: null, key: null, anchor: null };
  let hideTimer = null, lastPointer = "mouse";
  function renderTip() {
    const c = tipContent(tipState.key, tipState.target);
    if (!c) return false;
    const ref = REFS[c.ref] || {};
    const refHTML = ref.href
      ? `<a href="${esc(ref.href)}" data-role="ref">${esc(ref.label || ref.href)}</a>`
      : `<span data-role="ref">${esc(ref.label || "")}</span>`;
    tip.innerHTML = `<span class="cpuzzle-tt">${c.title}</span>${c.body}${c.value || ""}<span class="cpuzzle-tr">Defined on: ${refHTML}</span>`;
    tip.dataset.key = tipState.key;
    return true;
  }
  function place() {
    const rr = root.getBoundingClientRect();
    const k = (rr.width / (root.offsetWidth || rr.width)) || 1;        // reveal scale (transform or zoom)
    let a = tipState.anchor;
    if (!a) {
      const tr = tipState.target.getBoundingClientRect();
      a = { x: tr.left + tr.width / 2, top: tr.top, bottom: tr.bottom };
      if (tr.height > 4 * (G ? G.fs : 20) * k) a = { x: a.x, top: tr.top + 1.5 * (G ? G.fs : 20) * k, bottom: tr.top + 2.5 * (G ? G.fs : 20) * k };
    }
    tip.style.visibility = "hidden"; tip.hidden = false; tip.style.left = "0px"; tip.style.top = "0px";
    const tw = tip.offsetWidth, th = tip.offsetHeight, gapPx = 10;
    const cx = (a.x - rr.left) / k;
    const left = Math.max(0, Math.min(root.offsetWidth - tw, cx - tw / 2));
    const above = (a.top - rr.top) / k - th - gapPx;
    const below = (a.bottom - rr.top) / k + gapPx;
    // prefer above; go below when above would leave the widget's top by more than the bar height
    const top = above >= -bar.offsetHeight ? above : below;
    tip.style.left = left.toFixed(0) + "px"; tip.style.top = top.toFixed(0) + "px";
    tip.style.visibility = "";
  }
  function onDocKey(e) {
    if (e.key === "Escape" && tipState.open) {
      e.preventDefault(); e.stopPropagation(); if (e.stopImmediatePropagation) e.stopImmediatePropagation();
      const t = tipState.target; hideTip();
      if (t && root.contains(document.activeElement)) { /* keep focus where it was */ }
    }
  }
  function onDocDown(e) {
    if (!tipState.open || !tipState.pinned) return;
    if (tip.contains(e.target) || (tipState.target && tipState.target.contains(e.target))) return;
    hideTip();
  }
  function showTip(target, pinned, anchor) {
    if (!root.isConnected) return;
    clearTimeout(hideTimer);
    tipState.target = target; tipState.key = target.dataset.tip; tipState.pinned = !!pinned; tipState.anchor = anchor || null;
    if (!renderTip()) return;
    place();
    tipState.open = true;
    root.querySelectorAll("[aria-describedby]").forEach(n => n.removeAttribute("aria-describedby"));
    target.setAttribute("aria-describedby", tip.id);
    tipState.target.classList && tipState.target.classList.add("cpuzzle-tipped");
    window.addEventListener("keydown", onDocKey, true);
    document.addEventListener("pointerdown", onDocDown, true);
  }
  function refreshTip() { if (renderTip()) place(); }
  function hideTip() {
    clearTimeout(hideTimer);
    if (!tipState.open) return;
    tip.hidden = true; tipState.open = false; tipState.pinned = false;
    if (tipState.target) { tipState.target.removeAttribute("aria-describedby"); tipState.target.classList && tipState.target.classList.remove("cpuzzle-tipped"); }
    tipState.target = null;
    window.removeEventListener("keydown", onDocKey, true);
    document.removeEventListener("pointerdown", onDocDown, true);
  }
  const softHide = () => { if (!tipState.pinned) { clearTimeout(hideTimer); hideTimer = setTimeout(hideTip, 160); } };
  const tipTarget = n => (n && n.closest) ? n.closest("[data-tip]") : null;
  // Hover resolution for the dots. (1) Dots overlap where they are dense (the inverse-MSE row
  // near 1): of the dots under the pointer, the one whose centre is nearest wins, so every dot
  // can be reached by pointing at its centre. (2) The simple-average line and the median / p90
  // hit halos are drawn over the dots, and the replications near 1 (the ones whose verdict needs
  // full precision) sit under them; on those elements (class cpuzzle-hit) a dot under the
  // pointer wins, except on a median / p90 marker's visible bar.
  const isDot = n => !!(n && n.matches && n.matches('circle[data-tip="rel_rmse"]'));
  const isHit = n => !!(n && n.classList && n.classList.contains("cpuzzle-hit"));
  const nearestDot = e => {
    if (!document.elementsFromPoint) return null;
    let best = null, bd = Infinity;
    document.elementsFromPoint(e.clientX, e.clientY).forEach(n => {
      if (!isDot(n) || !root.contains(n)) return;
      const r = n.getBoundingClientRect(), d = Math.hypot(r.left + r.width / 2 - e.clientX, r.top + r.height / 2 - e.clientY);
      if (d < bd) { bd = d; best = n; }
    });
    return best;
  };
  const onBar = e => {
    const bar = e.target.parentNode && e.target.parentNode.querySelector("line:not(.cpuzzle-hit)");
    if (!bar) return false;                                      // simple-average line: defers everywhere
    const svg = bar.ownerSVGElement, k = svg.getBoundingClientRect().width / svg.viewBox.baseVal.width;
    return Math.abs(e.clientX - bar.getBoundingClientRect().left) <= 0.5 * (+bar.getAttribute("stroke-width")) * k + 1;
  };
  const pointerTarget = e => {
    if (isDot(e.target) || (isHit(e.target) && !onBar(e))) return nearestDot(e) || tipTarget(e.target);
    return tipTarget(e.target);
  };

  root.addEventListener("pointerdown", e => { lastPointer = e.pointerType || "mouse"; }, true);
  root.addEventListener("pointerover", e => {
    if (e.pointerType === "touch") return;
    if (tip.contains(e.target)) { clearTimeout(hideTimer); return; }
    const t = pointerTarget(e);
    if (!t || !root.contains(t)) return;
    if (tipState.pinned && tipState.target !== t) return;          // a pinned bubble stays until dismissed
    const anchor = t.tagName.toLowerCase() === "rect" ? { x: e.clientX, top: e.clientY - 6, bottom: e.clientY + 6 } : null;
    if (!(tipState.open && tipState.target === t)) showTip(t, false, anchor);
    else clearTimeout(hideTimer);
  });
  // moving within a dot or a hit element onto another dot fires no pointerover: re-resolve on move
  root.addEventListener("pointermove", e => {
    if (e.pointerType === "touch" || tipState.pinned) return;
    if (!isDot(e.target) && !isHit(e.target)) return;
    const t = pointerTarget(e);
    if (t && root.contains(t) && !(tipState.open && tipState.target === t)) showTip(t, false, null);
  });
  root.addEventListener("pointerout", e => {
    if (e.pointerType === "touch") return;
    const to = e.relatedTarget;
    if (to && (tip.contains(to) || (tipState.target && tipState.target.contains(to)))) return;
    softHide();
  });
  root.addEventListener("click", e => {
    if (tip.contains(e.target)) return;                             // links inside the bubble work normally
    const t = pointerTarget(e);
    if (!t || !root.contains(t)) return;
    const isBtn = t.classList.contains("cpuzzle-btn");
    if (isBtn && lastPointer !== "touch") return;                   // a mouse click on T_w just selects it
    if (tipState.open && tipState.pinned && tipState.target === t && !isBtn) { hideTip(); return; }
    const anchor = (t.tagName.toLowerCase() === "rect" && e.clientX) ? { x: e.clientX, top: e.clientY - 6, bottom: e.clientY + 6 } : null;
    showTip(t, true, anchor);
  });
  root.addEventListener("focusin", e => {
    const t = tipTarget(e.target);
    if (!t || !root.contains(t) || tip.contains(e.target)) return;
    if (t.classList.contains("cpuzzle-btn") && !(t.matches && t.matches(":focus-visible"))) return;
    showTip(t, false, null);
  });
  root.addEventListener("focusout", e => {
    const to = e.relatedTarget;
    if (to && tip.contains(to)) return;
    if (!tipState.pinned) hideTip();
  });
  tip.addEventListener("pointerenter", () => clearTimeout(hideTimer));
  tip.addEventListener("pointerleave", softHide);
  // Enter / Space on an info mark or legend item toggles a pinned bubble
  root.addEventListener("keydown", e => {
    const t = tipTarget(e.target);
    if (!t || t.classList.contains("cpuzzle-btn")) return;
    if (e.key === "Enter" || e.key === " ") {
      e.preventDefault();
      if (tipState.open && tipState.pinned && tipState.target === t) hideTip(); else showTip(t, true, null);
    }
  });
  // window.Reveal, never the bare identifier: OJS resolves a free name as a cell / global and
  // throws when it is absent (the notes have no Reveal).
  const RV = typeof window !== "undefined" ? window.Reveal : null;
  if (RV && RV.on) RV.on("slidechanged", () => { if (tipState.open) hideTip(); });

  // ---- keyboard: radio-group arrows; keep reveal.js from also turning the slide ----
  // PageUp / PageDown are deliberately NOT stopped: there is no slider here, and a presentation
  // clicker sends them; focus stays on a T_w button after a click, so stopping them would
  // freeze the deck. Escape is handled by the bubble (only while one is open).
  const NAV = ["ArrowLeft", "ArrowRight", "ArrowUp", "ArrowDown", "Home", "End", " ", "Enter"];
  grp.addEventListener("keydown", e => {
    let k = null;
    if (e.key === "ArrowRight" || e.key === "ArrowDown") k = Math.min(TW.length - 1, j + 1);
    else if (e.key === "ArrowLeft" || e.key === "ArrowUp") k = Math.max(0, j - 1);
    else if (e.key === "Home") k = 0;
    else if (e.key === "End") k = TW.length - 1;
    if (k !== null) { e.preventDefault(); e.stopPropagation(); setJ(k); BTN[k].focus({ preventScroll: true }); if (tipState.open) showTip(BTN[k], false, null); }
  });
  root.addEventListener("keydown", e => { if (NAV.indexOf(e.key) >= 0) e.stopPropagation(); });
  root.addEventListener("keyup", e => { if (NAV.indexOf(e.key) >= 0) e.stopPropagation(); });

  // ---- first paint and re-layout on resize ---------------------------------------------
  build(); buttons(); words(); draw(mkCur);
  if (typeof ResizeObserver !== "undefined") {
    const ro = new ResizeObserver(() => {
      if (!root.isConnected) { ro.disconnect(); hideTip(); return; }
      const w = root.clientWidth, host = parseFloat(getComputedStyle(root).fontSize) || 16;
      if (w > 0 && (Math.abs(w - (G.lastW || 0)) > 2 || host !== G.lastHost || inDeck() !== G.deck)) {
        const was = tipState.open; hideTip(); build(); draw(mkCur); words();
        if (was) { /* targets were rebuilt; the bubble closes rather than float over a stale element */ }
      }
    });
    ro.observe(root);
  }

  // hooks for headless QA (not part of the API)
  root._cpuzzle = {
    get Tw() { return TW[j]; }, get shown() { return TW[jShown]; }, get geom() { return G; },
    get animating() { return !!raf; }, get tip() { return { open: tipState.open, pinned: tipState.pinned, key: tipState.key }; },
    refs: REFS, set: k => setJ(TW.indexOf(k)), reset: () => setJ(j0), hideTip, data: D
  };
  return root;
}