Conformal intervals become extreme/null after inverse log transform (expm1) — what's the recommended approach?

[gbgoutha]

What happened + What you expected to happen

Conformal intervals produce extreme/null values after inverse log transform (expm1)

Description

When using ConformalIntervals with AutoARIMA on log1p-transformed data, the conformal scores (computed in log-space) produce extreme or null hi values after applying expm1 to inverse-transform back to original scale.

This is because conformal scores are additive in log-space but become multiplicative in original space after expm1.

Example

column	log-space value	after expm1
point	6.03	414
lo-70	5.31	202
hi-70	11.88	144,764

The hi is 350x the point forecast — clearly nonsensical for a demand prediction interval.

Root cause

In _conformity_scores (models.py L122-148), the absolute residuals |mean - y_test| are computed in the transformed (log) space. In _add_conformal_distribution_intervals (models.py L69-96), these scores are added/subtracted from the point forecast, still in log-space:

scores = np.vstack([mean - cs, mean + cs])  # all in log-space

### Versions / Dependencies

<details><summary>Click to expand</summary>
Dependencies:
statsforecast 2.0.3
</details>


### Reproducible example

```python
from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA
from statsforecast.utils import ConformalIntervals
import numpy as np
import pandas as pd

np.random.seed(42)
n = 260  # 5 years weekly
y = np.random.lognormal(mean=6, sigma=0.5, size=n)
df = pd.DataFrame({
    'unique_id': ['series_1'] * n,
    'ds': pd.date_range('2020-01-06', periods=n, freq='W-MON'),
    'y': y
})

# Log transform before fitting
df['y'] = np.log1p(df['y'])

sf = StatsForecast(models=[AutoARIMA(season_length=52, D=1)], freq='W-MON')
ci = ConformalIntervals(h=26, n_windows=3)

fcst = sf.forecast(df=df, h=26, prediction_intervals=ci, level=[70, 90])

# Inverse transform — hi values explode
for col in ['AutoARIMA', 'AutoARIMA-lo-70', 'AutoARIMA-hi-70',
            'AutoARIMA-lo-90', 'AutoARIMA-hi-90']:
    fcst[col] = np.expm1(fcst[col])

print(fcst[['AutoARIMA', 'AutoARIMA-lo-70', 'AutoARIMA-hi-70']].describe())
# hi-70 mean will be orders of magnitude larger than the point forecast

Issue Severity

High: It blocks me from completing my task.

[eskayML]

looks like the conformal scores are being computed in log-space but then transformed back to original scale using expm1, which is causing the extreme/null values. one possible fix could be to compute the conformal scores in the original scale from the start, or to apply a log transform to the scores before adding/subtracting them from the point forecast. should i open a PR for this?

[gbgoutha]

Thanks for the suggestion! Yes, that's exactly the root cause — the conformal scores are computed in log-space (since the transformation is applied externally before fitting AutoARIMA), but the back-transformation via expm1 is applied incorrectly, leading to extreme or null values in the prediction intervals.

The blambda parameter inside AutoARIMA would ideally handle this internally, but it currently only works when set to None. Setting blambda=0 (log transform) is broken, so the transform has to be applied outside the model for now.

A PR to fix this would be welcome. The cleanest approach would likely be to ensure conformal scores are either:

• Computed entirely in the original scale (back-transform predictions before scoring), or
• Kept in log-space throughout and back-transformed only after the full interval is constructed

The latter is generally more numerically stable and preserves the asymmetry of the log-normal distribution correctly. Happy to review a PR if you'd like to go ahead!

[Whatsonyourmind]

The root cause is well-identified. To add some mathematical framing: conformal scores need to be equivariant under the transform. For a monotone transform g(x) = log1p(x), the correct approach is to construct the prediction interval in log-space and then transform the entire interval, not transform the scores additively.

Concretely, the fix should change the interval construction from:

# Current (broken): additive scores in log-space
lo = mean_log - score
hi = mean_log + score
# Then expm1(lo), expm1(hi) -- explodes because expm1 is convex

to:

# Correct: transform the full interval endpoints
lo = expm1(mean_log - score)
hi = expm1(mean_log + score)

But even this will produce asymmetric intervals in the original scale (wider on the high side), which is actually correct for lognormal-like data. The issue reporter's 350x ratio is extreme, but some asymmetry is expected since demand data is right-skewed.

The numerically stable approach for the PR: apply the conformal procedure entirely in log-space, construct the quantile distribution in log-space, and let the user apply expm1 after. The key principle is that conformal coverage guarantees hold regardless of the scale, so the intervals in log-space have valid coverage -- you just need to back-transform the endpoints (not the scores) to get valid intervals in the original scale.

One thing to watch for: if blambda=0 is truly broken (as @gbgoutha suggests), fixing conformal intervals for external transforms might be a band-aid. The cleaner long-term fix is making blambda=0 work so that the model handles the transform internally and conformal scores are computed on the right scale automatically.