# Assignment 4: An h-step forecaster for a plant channel
**Module:** Week 04 · **Released:** L8 (2026-09-21) · **Due:** 2026-09-28 · **6 points on Canvas.**

## Goal
Forecast the stripper temperature of the Tennessee Eastman plant at five horizons, from 3 minutes to 2 hours, and show honestly whether your model is worth using. "Worth using" has a precise meaning here: better than the two free forecasts, persistence and the mean, on simulation runs the model never saw.

Everything you need is in [Lecture 8](../../lectures/l08/notes.md). The demo notebook does the same steps on reactor pressure; this assignment does them on a different channel, with a validation step the demo skips.

**You should be able to:**
- build a forecasting table whose features are known at time `t` and whose target is `y[t+h]`, inside each run;
- score persistence and the mean at each horizon, and compute a skill score against the better one;
- fit a direct forecaster per horizon and a recursive one, and compare them;
- choose hyperparameters on validation runs, never on test runs;
- show how much a shuffled split flatters a model on a single series, using `TimeSeriesSplit` with a gap.

## The data
The fault-free training file from [Rieth et al. (2017)](https://doi.org/10.7910/DVN/6C3JR1), the same simulations Lecture 8 uses. The course hosts it as Parquet, 25 MB, so you need no R reader:

```bash
mkdir -p data && cd data
curl -fLO https://kitchin-services.cheme.cmu.edu/f26-06763/data/tep_fault_free_training.parquet
curl -fLO https://kitchin-services.cheme.cmu.edu/f26-06763/data/SHA256SUMS
shasum -a 256 --ignore-missing -c SHA256SUMS
```

`pl.read_parquet` opens it directly. The original is the `TEP_FaultFree_Training.RData` file in the Dataverse record above, converted with no other change. **Do not commit the data.**

| | |
|---|---|
| Runs | 500 independent simulations, `simulationRun` 1 to 500 |
| Samples | 500 per run, `sample` 1 to 500, every 3 minutes (25 hours) |
| Target | `xmeas_18`, stripper temperature, °C |
| Inputs you may use | the target's own past, and the 11 valve positions `xmv_1` to `xmv_11` **at time `t` only** |

## Fixed choices
These are fixed so that every submission can be checked against the same numbers.

| Choice | Value |
|---|---|
| Horizons `h` | 1, 5, 10, 20, 40 samples (3, 15, 30, 60, 120 minutes) |
| Training runs | 1 to 300 |
| Validation runs | 301 to 400, for choosing hyperparameters only |
| Test runs | 401 to 500, touched once, at the end |
| Forecast origins | every `t` from 40 to `500 - h` in each test run |
| Mean baseline | the mean of `xmeas_18` over **all** samples of runs 1 to 300 |
| Persistence | `y[t+h] = y[t]` |
| Error | RMSE in °C |

Starting every origin at `t = 40` means every model and every baseline is scored on exactly the same rows, whatever lag depth you pick, as long as it is 40 or fewer.

## Tasks

### 1. Build the forecasting table
Write one function that returns, for a given `h` and lag depth, one row per `(run, t)` with the features known at `t` and the target `y[t+h]`. Every shift must run inside a run (`.over("simulationRun")` in Polars, or a per-run loop). Drop rows whose lags or target fall outside their run. No feature may come from after `t`.

### 2. Score the baselines
On the test origins, compute persistence and mean RMSE at every horizon.

### 3. Choose hyperparameters on the validation runs
Fit a scaled ridge regression (`make_pipeline(StandardScaler(), Ridge(alpha=...))`) on the training runs at `h = 10`, and score it on the validation runs for at least three settings: vary the lag depth, and the penalty if you like. Write every setting you tried to `results/validation.csv`. Use the best setting for everything that follows. **The test runs play no part in this choice.**

### 4. Fit direct and recursive forecasters
- **Direct:** one model per horizon, trained on runs 1 to 300 with target `y[t+h]`. Lags of `xmeas_18` plus the valves at `t`.
- **Recursive:** one model for `y[t+1]` from lags of `xmeas_18` only, applied `h` times, feeding each prediction back into the lags. Say in your report why this one does not use the valves.

You may add a second direct model of your own choice (gradient boosting, a random forest, extra window features), as long as the ridge results are also reported.

### 5. Save predictions and the results table
Write `results/predictions.parquet` with one row per test origin and horizon:

| Column | Type | Meaning |
|---|---|---|
| `run` | int | `simulationRun`, 401 to 500 |
| `t` | int | the forecast origin, `sample` at time `t` |
| `h` | int | horizon in samples |
| `y_true` | float | `xmeas_18` at sample `t + h` of that run |
| `persistence` | float | your persistence forecast |
| `mean` | float | your mean forecast |
| `direct` | float | your direct forecast |
| `recursive` | float | your recursive forecast |

Write `results/forecast.csv` with one row per horizon and at least these columns: `h`, `persistence_rmse`, `mean_rmse`, `direct_rmse`, `recursive_rmse`, `skill`, where `skill` is `1 - direct_rmse / min(persistence_rmse, mean_rmse)`.

### 6. Shuffle, then don't
Take runs 1 to 10, **one run at a time**, at `h = 10`. For each run, score a random forest (or another flexible model) two ways, with the same features as your direct model:

- `KFold(n_splits=5, shuffle=True, random_state=0)`
- `TimeSeriesSplit(n_splits=5, gap=10)`

Score persistence on the same test folds. Write `results/leaky.csv` with columns `run`, `split` (`shuffled` or `time`), `h`, `model_rmse`, `persistence_rmse`, one row per run and split.

### 7. Write `REPORT.md`
**Two pages maximum**, six sections in this order:

1. **The table.** What a row holds, what is known at `t`, and how you kept each shift inside its run.
2. **Baselines and skill.** Your `forecast.csv` as a table. Where persistence stops beating the mean, and whether that matches the autocorrelation of `xmeas_18`.
3. **Validation.** What you tried, what you chose, and how much it mattered.
4. **Direct against recursive.** Which held up better, where, and why. Why the recursive model has no valves.
5. **Shuffle against time.** Your mean `leaky.csv` numbers, which one you would report to a plant manager, and why the gap appears.
6. **Limits.** One situation where this forecaster would give a confident wrong answer in a real plant.

End with one line disclosing generative-AI use.

## Names to use

| Thing | Name |
|---|---|
| Code | any `.py` files or notebooks in the project |
| Data | `data/tep_fault_free_training.parquet`, not committed |
| Validation sweep | `results/validation.csv` |
| Predictions | `results/predictions.parquet` |
| Results table | `results/forecast.csv` |
| Shuffle comparison | `results/leaky.csv` |
| Report | `REPORT.md` |

## Submit
From your project root, after running your code:

```bash
uv run --no-project https://kitchingroup.cheme.cmu.edu/f26-06763/a04-evidence.py \
    --andrew-id yourid --name "Your Name"
```

Upload the resulting `evidence.pdf` to Canvas.

- **Keep `--no-project`**, so a broken `uv.lock` in your project cannot stop the script.
- The script reads the Rieth file from your `data/` directory (pass `--data` if it is elsewhere) and **recomputes** the baselines and the targets itself. It does not re-run your code, and it does not download anything.
- **Read the PDF before uploading.** It prints your score by group. If a check fails, fix it and rerun, and submit anyway if you run out of time.
- The PDF prints the script's sha256, which matches <https://kitchingroup.cheme.cmu.edu/f26-06763/a04-evidence.py.sha256> when run from the URL.

## Grading
Five groups are scored by the script and the sixth by your TA. Within a group, checks are equally weighted.

| # | Group | Pts | Checks |
|---|---|---|---|
| 1 | **The table** | 1.0 | `predictions.parquet` with the eight columns; every test run and origin present for every horizon, and no other runs; `y_true` is really `xmeas_18` at `t + h`; shifts kept inside each run in your code |
| 2 | **Baselines** | 1.0 | `forecast.csv` has the five horizons; persistence and mean match the script's recomputation; `skill` matches its definition |
| 3 | **Models** | 1.25 | direct and recursive RMSE recomputed from your predictions match `forecast.csv`; direct beats the better baseline at four or more horizons; recursive differs from direct beyond `h = 1`; a scaler inside a pipeline in your code |
| 4 | **Honest evaluation** | 1.25 | `validation.csv` with three or more settings; validation uses runs 301 to 400 in your code; `leaky.csv` with both splits for runs 1 to 10; `TimeSeriesSplit` with a gap of at least 10 in your code; shuffled `KFold` in your code |
| 5 | **Report numbers** | 0.5 | the RMSE values quoted in `REPORT.md` appear in `forecast.csv`; the report quotes both `leaky.csv` means |
| 6 | **REPORT.md** | 1.0 | the six sections and the AI-use line, read by your TA |

A check the script cannot decide is held for your TA rather than lost. Your TA can adjust in either direction when the work is clearly better or worse than the checks show.

## AI use
Generative AI is allowed with disclosure in `REPORT.md`. Editing the generated PDF by hand is falsifying a submission. You must be able to explain your table, why the gap has to be at least `h`, and why your shuffled score is lower.

## Stretch (not graded)
- **A residual detector.** Download the faulty training file (datafile `3031242`, about 500 MB). Set a threshold on your one-step residuals from validation runs, then run it on fault 1 and one other fault of your choice. Report the false-alarm rate on the test runs and the detection delay on each fault. Faults are introduced one hour into each faulty training run.
- Plot the autocorrelation of `xmeas_18` and predict the persistence-mean crossover from it before looking at your table.
- Try `skforecast`'s `ForecasterDirect` and `ForecasterRecursive` and check they reproduce your numbers.