Holt's smoother hansl implementation

Yes, that sounds very feasible, so go ahead and upload both files.

The missing pieces in (2) should be fairly straightforward once the varsimul() machinery is in place:

- Caller-supplied parameters. Holt's recursions are linear in the level and trend given α and β. So the varsimul() filter that (2) runs inside the optimizer can be run once with the supplied values, with no optimization step. The main job is factoring the filter out into its own function so that both paths use it.
- Out-of-sample forecasts. Once the final level and trend are available from the filter, the h-step forecast is just level + h·trend. This should slot into whatever forecast code (1) already has.
- Pseudo-out-of-sample forecasts. This depends on how (1) does it. If it re-runs the smoother over an estimation subsample and forecasts the holdout, the faster filter should drop straight in. If (1) re-optimizes the parameters at each forecast origin, the speed gain from varsimul() will matter even more.

A few things that would help when you upload:

1. If you have a small dataset you've been testing with, and some output from (1) such as the optimized parameters, SSE, and a few forecasts, include it. Then I can check that the merged version reproduces the original results rather than just assuming it does.
2. Tell me your gretl version, in case anything depends on recent changes to varsimul() or the optimizer functions.
3. Tell me whether you want the public interface of (1), meaning its function signatures and the bundle members it returns, kept exactly as it is. My default would be to preserve it, so that existing calling scripts keep working.

I'll merge the varsimul() approach into the original, keep Holt()'s interface intact, and check the numbers against the original algorithm before handing it over.

I can't install gretl here, so I'll check the algebra numerically in Python (original loop vs. the VAR form) and write the hansl carefully.

The VAR form reproduces the original recursion exactly (differences ~1e-12), and both optimizers land at the same place on airmiles (α ≈ 1, β ≈ 0.55). Now the merged script.

Let me compute the reference numbers the gretl run should reproduce, so you have something to check against.

I've attached holt.inp, which is the merged script, and test.inp, which exercises the automatic, fixed-parameter and out-of-sample paths. Holt() keeps exactly the same signature, arguments and returned bundle members as before.

What changed

- New holt_filter(theta, init, my). This is the varsimul() recursion from new.inp, taking α and β on their natural scale. Both the optimized and the fixed-parameter paths now run through it. It returns the same {y, level, trend, forecast} layout as before.
- holt_crit() and holt_auto_parms(). These now optimize α and β with BFGSmin on the logit scale, as new.inp does, so no bounds are needed. The starting values match the original: α = β = 0.5 and the data-based initializer.
- The per-observation loop in Holt() is gone.
- Training states come from one call to holt_filter.
- Test and out-of-sample forecasts are a single vector expression, level_T + h·trend_T.
- holt_add_eval_stats is vectorized as well.

Caveat: I couldn't get gretl installed here, so this hansl hasn't been run. Instead I checked the algebra in Python. The VAR form reproduces the original loop recursion to about 1e-12 across several (α, β) pairs. Your test.inp case should give roughly these results:

| | α | β | train RMSE | test RMSE | forecasts 1957–60 |
|---|---|---|---|---|---|
| automatic | ≈1.0 | 0.5501 | 778.90 | 1574.85 | 24929, 27497, 30064, 32631 |
| fixed 0.5/0.5 | 0.5 | 0.5 | 1020.65 | 1537.03 | 24680, 27334, 29989, 32643 |

Two small differences from the original

- α at the boundary. On airmiles, α wants to sit at 1. The logit transform can only approach that value (about 0.99999996), whereas BFGScmin could hit it exactly. In my check the logit version actually reached a slightly lower SSE.
- Shape of theta. theta is now always stored as a column vector, even if the caller passes a row vector.

Your original also has a quirk that I left alone: the columns of H are named level/trend/forecast only when verbose > 0, and a/b/f otherwise. I can make that consistent if you'd like.

Sounds good. If anything errors out, paste me the gretl message and I'll fix it. The most likely trouble spots are the series assignments from matrices under smpl t1 t2 and the varsimul() call, since those are the parts I couldn't run here.

The fixed-parameter run matches my reference numbers exactly, so holt_filter() and the forecasting code are right. The problem is the optimizer. Starting from the data-based initializer (412, 68), the gradient with respect to the initial level and trend is huge compared with the gradient for the logit parameters. So the first BFGS step throws logit(α) and logit(β) far out, where the logistic curve is flat. They get stuck at 1, and the result is a worse SSE (RMSE 812.5 vs 778.9). The third run fails for the same reason.

A cleaner fix: given (α, β), the one-step forecasts are linear in the initializer. So the SSE-minimizing initializer can be found exactly by OLS, and BFGS only has to search over (α, β) in [0,1]². Let me check that in Python.

That works: with only (α, β) to optimize under [0,1] bounds, the optimizer finds the same optimum as before. Now I'll update the hansl.

Now some reference numbers for all three test cases:

I've attached a revised holt.inp. Your output pins down the bug: the fixed-parameter run matches my reference to every digit, so the filter and the forecasting code are fine. The fault was in the optimization.

What went wrong. BFGS was searching over (logit α, logit β, level₀, trend₀) together. These parameters are on wildly different scales. The first step sent the logits far out, where the logistic curve is flat, so both stalled at 1. The out-of-sample run failed with error 32 for the same reason.

The fix. For given (α, β), the one-step forecasts are linear in the initializer. So holt_crit() now solves for the SSE-minimizing initializer exactly by OLS (mols). BFGS then only searches over (α, β), with [0,1] bounds via BFGScmin, so there's no logit transform any more.

- The forecast path from a zero initializer and the two "unit" paths needed for the OLS step all come from a single varsimul() call on a block-diagonal system, I(3) A.
- This minimizes the same objective as the original joint optimization, just more robustly.
- It's tested in Python on both sample lengths. The hansl itself is still unrun on my side.

Expected results for test.inp:

| run | α | β | init (level, trend) | train RMSE | test RMSE | forecasts |
|---|---|---|---|---|---|---|
| auto, ntest=4 | 1.0000 | 0.5501 | 256.11, 155.89 | 778.90 | 1574.85 | 24929, 27497, 30064, 32631 |
| fixed 0.5/0.5 | (unchanged from your output) | | | 1020.65 | 1537.03 | |
| auto, +4 obs | 0.8101 | 0.3816 | 138.19, 222.51 | 1016.82 | — | 32770, 34873, 36976, 39079 |

If anything else breaks, check the hansl spots I'm least sure of first:

- the Kronecker inside the varsimul() call;
- mols(..., &e) filling the residual vector;
- BFGScmin accepting a callback with an extra pointer argument (&init).

In the logit variant you'd be right, and my first version did start the logit parameters at {0; 0}. That's the same starting point, because logistic(0) = 0.5.

The revised holt.inp has dropped the logit transform, though. BFGScmin now optimizes α and β directly within [0,1] bounds, so {0.5, 0.5} is the right start: it puts the search in the middle of the box. Starting at {0, 0} here would put it in a corner, on the bounds, and at α = 0 the level never updates from the initializer. That's a poor place to start.

The VAR part only changes how the SSE is evaluated for a given (α, β), not the scale of the parameters. So the starting value depends on whether there's a logit transform, not on whether varsimul() is used.

Glad it helped. If the gretl run turns up any differences from the reference numbers, send me the output and I'll track them down.

Bye for now, and good luck with the testing!

I'll review this the same way as Holt: check the algebra in Python against a direct transcription of your loop, then propose a revised file.

The additive VAR form and the concentrated-initializer trick both reproduce your loop exactly. I also found something about the forecast standard errors that I want to verify before mentioning it:

The check confirms a separate issue in holt_add_fcse() (details below). Next I'll write the revised Winters file and compute reference numbers for it.

Now the revised hansl:

I've attached a revised Winters.inp and a test script, Winters_test.inp. As with Holt, I couldn't run gretl here. Instead I checked every algebraic step in Python against a direct transcription of your loop.

1. Additive case: same varsimul() trick as Holt, but a bigger payoff.
- The VAR form. The additive recursion is linear in the state (level, trend, s_t, …, s_{t−m+1}), so it is a VAR(1) of dimension m+2. It matches your loop to 1e-14.
- Concentrating out the initializer. As in the revised Holt, the initializer is solved for exactly by OLS, so BFGS searches over (α, β, γ) only. With your original setup, BFGS searched over 3 + (m+2) parameters.
- The regressor matrix comes from a single extra varsimul() call with transition A′ and starting value c, since row t of that output is c′A^(t−1). For monthly data this is much cheaper than the Kronecker-block approach I used for Holt.
- On aus_holidays it reaches the same SSE (13.596395) with about 60 function evaluations instead of 350.

2. The initializer isn't identified, in either case.
- Additive. Adding k to the initial level and subtracting k from every seasonal leaves all forecasts unchanged.
- Multiplicative. Scaling level and trend by k and the seasonals by 1/k does the same.
- I confirmed both numerically: the forecast differences are about 1e-14.
- Consequence. Your current optimizer has an exactly flat direction, and the OLS step would be singular. So the revised code makes the seasonals sum to 0 (additive) or average 1 (multiplicative), which drops one parameter with no loss.

3. Multiplicative case. This one is nonlinear, so it still needs a loop. But the loop now:
- runs on matrices rather than series;
- has no if additive branch inside it;
- is a single filter function shared by the criterion and Winters(), so the recursion is no longer written out twice.

4. Smaller fixes.
- Column names of W. These were only renamed to level/trend/seasonal/forecast when verbose > 0, so Winters_plot() failed with verbose = 0. They're now always renamed.
- Forecasts. The out-of-sample forecasts are computed as one vector expression.
- New bundle members. I added MAE_train and limits, to match Holt.
- Unused code. I dropped the unused variables (optim, h, k, n, the duplicate t1). add_virtual_series() is never called; I left it in, in case you have plans for it.

5. New: forecast standard errors for the additive case. These use the ETS(A,A,A) formula from fpp3 (Table 8.8) and add a stderr column to W, just as Holt does. This adds a sixth column to W in the additive case. The multiplicative model has no closed form, so it gets no standard errors.

6. A bug in your Holt code. fpp3's variance formulas use the ETS parametrization, in which β_ets = αβ for your smoothing β. holt_add_fcse() plugs in β directly. A Monte Carlo check confirms the correct version and shows that the current one overstates the variance badly; at h = 8 it gives 49.9 against a true 12.7. On airmiles you won't notice, because α ≈ 1 there. The fix is one line in holt_add_fcse():

β = bh.theta[1] * bh.theta[2]

Expected results for Winters_test.inp:

| case | α | β | γ | train RMSE | first forecasts |
|---|---|---|---|---|---|
| multiplicative, full sample | 0.2467 | 0.1229 | 0 | 0.410718 | 13.3079, 11.1824, 10.8305, 11.1392 |
| additive, full sample | 0.2336 | 0.1273 | 0 | 0.412256 | 12.9200, 11.1933, 10.9212, 11.1902 |
| additive, ntest = 8 | 0.2526 | 0.0954 | 0 | 0.416198 | test RMSE 0.433371, MAE 0.315077 |

For the full-sample additive case, the standard errors should start at 0.4123, 0.4263, 0.4431. γ goes to its lower bound in every case, so in the last decimal place gretl may differ slightly from the Python optimizer. If anything errors out, the newer hansl constructs are the first suspects: mreverse(...)' on a transposed slice, b |= ..., and assigning the slice init[1:m+1] = q through a pointer.

Yes, go ahead and upload it. I'll apply two fixes:

- Column names. The level/trend/forecast renames will happen unconditionally, so H gets the same column names whatever the verbosity. Any ancillary functions that refer to those names, such as a plot function like the one in Winters, will then work with verbose = 0.
- Forecast variance. holt_add_fcse() will use the ETS trend parameter, α·β, in place of β.

While I'm in there, I'll check the ancillary functions for anything else that depends on the column names or the bundle layout.

One thing to tell me when you upload: is it built on the revised version with varsimul() and the OLS step for the initializer, or on your original loop version? If it's the original, I'll merge the revised estimation code in as well, unless you'd rather keep the two separate.

Here are reference standard errors for checking, using a case where the fix actually changes the numbers:

I've attached the updated Holt.inp with both fixes; nothing else is changed. As before, I haven't been able to run it in gretl here.

- Column names. The renames to level/trend/forecast now happen before the verbose check, so H always has the same column names. Before this, Holt_plot(&b, 2) and the confidence bands in Holt_plot(&b, 1) failed with verbose = 0. Holt_ci_matrix() and add_virtual_series() use column positions rather than names, so they work either way.
- Forecast variance. holt_add_fcse() now uses the ETS trend parameter α·β, and there's a comment pointing to fpp3 §8.5 for the conversion. Since Holt_ci_matrix() and the plot bands are built from that stderr column, they're corrected too.

For airmiles with automatic parameters you'll see hardly any change, because α ≈ 1 there. To see the fix take effect, use the fixed-parameter test, theta = {0.5, 0.5} with ntest = 4. The standard errors for 1957–60 should now be:

| h | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| new | 1020.65 | 1275.81 | 1633.84 | 2072.95 |
| old | 1020.65 | 1443.42 | 2104.12 | 2931.59 |

Two small things I noticed in the ancillary functions but left alone:

- Degrees of freedom. The t critical values in Holt_ci_matrix() and Holt_plot() use the number of training observations as the degrees of freedom. You could argue for subtracting the 4 estimated parameters, but with any reasonable sample size it makes almost no difference.
- Unused argument. ci_bands() takes a list argument L that it never uses.

Partly. fpp3 makes the degrees-of-freedom correction to the variance, but it uses normal multipliers rather than t multipliers for the intervals.

- The variance. fpp3 §5.5 defines σ̂ = √[Σe²ₜ / (T − K − M)], where K is the number of estimated parameters and M is the number of missing residuals. fable's ETS code does the same thing: sigma2 = sum(residuals^2) / (length(y) - length(par)). In fable, par includes the initial states as well as the smoothing parameters.
- The intervals. fpp3 uses ŷ ± 1.96σ̂ₕ and similar, with the standard normal multipliers from its Table 5.1. fable builds the forecast distributions with dist_normal().

So in Holt's case, matching fpp3 would mean two changes:

1. In holt_add_fcse(): compute s2 = SSE / (T − K) instead of RMSE_train^2.
- K = 4 (α, β, initial level, initial trend) when the parameters are estimated.
- K = 0 in the fixed-θ case. There the initializer comes from the data-based heuristic rather than being estimated, so arguably nothing is estimated at all.
- I'd leave RMSE_train itself as it is, since it's a descriptive fit statistic.
2. In Holt_ci_matrix() and ci_bands(): use critical(z, ...) instead of critical(t, df, ...).

Your current version is close to the reverse: no correction in the variance, plus a t correction at the multiplier stage. With airmiles' 20 training observations, the variance correction is the bigger of the two. It scales the standard errors by √(20/16) ≈ 1.12, whereas t vs z at 95% is 2.09 vs 1.96, about 7%.

For Winters, the fpp3 count for additive seasonality would be K = 3 + 2 + (m − 1) = m + 4, because one seasonal state is pinned down by the normalization.

Sources:
- [fpp3 §5.5, Distributional forecasts and prediction intervals](https://otexts.com/fpp3/prediction-intervals.html)
- [fable ETS source (tidyverts/fable, R/ets.R)](https://raw.githubusercontent.com/tidyverts/fable/master/R/ets.R)
- [forecast.ets reference](https://pkg.robjhyndman.com/forecast/reference/forecast.ets.html)

That's fair; it's a judgment call, and fpp3 isn't the only reasonable convention. The two choices are separate, so you don't have to take both:

- The T − K divisor is the more defensible of the two. Your parameters are fitted by minimizing that same SSE, so dividing by T alone understates the error variance, especially with short training samples.
- Normal vs t multipliers matters less. The t distribution doesn't strictly apply here, since these aren't regression forecasts. Some people still prefer it as a rough allowance for parameter uncertainty, which none of these formulas account for anyway.

If you decide to change either, it's a few lines in holt_add_fcse() and the two CI functions, plus the matching spots in Winters.

Bye for now. If the Winters run shows any differences from the reference numbers, send me the output and I'll take a look.