Probabilistic solar forecasting: confidence interval, prediction interval and PICP

A solar irradiance forecast needs a margin. At Ajaccio, in 2007, one hour ahead, the mean value forecast by the autoregression is known to ±0.5 W/m² (confidence interval), but the margin that aims to contain 95% of the next measurements is ±174 W/m² (prediction interval; 95.2% measured in 2007). This page shows why the two margins differ, why “forecast ± 2 RMSE” goes wrong, and how conformal prediction sets the margin from data.

An experiment from the solar forecasting course. The numbers are recomputed here from the measurements; they are not results from a paper. Interval method: Voyant et al., Solar Energy 2025 and JRSE 2022. French version.

Try it

Each black dot is an hourly measurement of GHI (global horizontal irradiance, the solar radiation received by a horizontal surface, in W/m²). The blue line is the forecast made one hour before. The light blue band is the prediction interval: where the next measurement should fall. The thin orange band, often invisible because it is so narrow, is the confidence interval of the mean forecast. Red dots are measurements that fell outside the band.

Model? Prediction interval? Nominal level? 95 % Date (2007) Window? Confidence interval?

measurement · measurement outside the prediction interval · 1 h forecast · prediction interval · confidence interval

PICP, window?

PICP, whole of 2007?

Mean band width, 2007?

CI: mean half-width?

CRPS, 2007?

PICP by hour of the day, 2007 (end hour of the measurement). Orange line: nominal level.

PICP by type of day, 2007, and mean band width (W/m²) under each bar.

Standard deviation of the series or RMSE of the model
Mean daytime GHI, 2007?
Standard deviation of the series $σ_{G H I}$ , 2007?
Model RMSE, calibration 2006?
Model RMSE, test 2007
Half-width of the chosen interval?

Two intervals, two questions

Confidence interval (CI): It answers “where is the mean of the measurements we would get in this situation?”. It only holds the uncertainty on the model parameters, estimated from a finite amount of data. It shrinks as data are added: with 8 years of training, it is ±0.5 W/m² here for AR and ±1.9 W/m² for ELM. Persistence has no estimated parameter: its CI reduces to the line.
Prediction interval (PI): It answers “where will the next measurement fall?”. It adds to the CI the spread of the measurements around their mean: a cloud passes or not. This part does not vanish with more data. It is the interval that matters to run a grid or a battery.
Nominal level 1 − α: The share of measurements the prediction interval aims to contain, for example 95% (α = 0.05, the share left outside).

The same distinction on a simple regression line: confidence or prediction interval?

Standard deviation of the series or error of the model?

Two numbers with similar names that have nothing in common:

Standard deviation of the series, $σ_{G H I}$: The spread of the measurements around their mean. Daytime 2007: mean 372 W/m², $σ_{G H I}$ = 282 W/m². It mixes nightfall, seasons and clouds: everything a model already knows how to forecast.
Residual e: Measurement minus forecast: $e = G H I - \hat{G H I}$ , in W/m².
RMSE: Root mean square error: the square root of the mean squared residual. It is the spread of the measurements around the forecast. If the residuals average zero (no bias), the RMSE is exactly their standard deviation; otherwise, ${R M S E}^{2} = {standard deviation}^{2} + {bias}^{2}$ . For AR, on the 2006 calibration year: bias -6.6 W/m², residual standard deviation 77.4 W/m², RMSE 77.7 W/m². The bias is small, so the two almost coincide.

Forecast ± 2 $σ_{G H I}$ would give a band of ±564 W/m²: almost the whole sky. Forecast ± 2 RMSE (exactly 1.96 RMSE at 95%) is the classical recipe. It assumes two things: Gaussian residuals (bell-shaped, without heavy tails) and homoscedastic residuals (the same spread at every hour and in every weather). Neither holds for solar irradiance:

over the whole of 2007, the AR ±1.96 RMSE band contains 93.4% of the measurements instead of 95%;
by hour, from 86% (hour 12, around noon, when irradiance is high and clouds are costly) to 98%; and at least 98% at sunrise and sunset hours, where errors are small;
by type of day: 98.4% on clear days (the band is too wide), 85.2% on variable days (too narrow), 89.7% on overcast days.

Select “Gauss ±z·RMSE” above and look at the two small charts: the bars move away from the target line. The residual histogram also shows a peak too sharp and tails too heavy for a Gaussian bell. To explore departures from normality in skewness and tail weight, try the Jarque–Bera test.

Formulas: standard deviation, RMSE and Gaussian interval

$σ_{G H I}^{2} = \frac{1}{N} \sum_{t = 1}^{N} {({G H I}_{t} - G H I)}^{2}$

${R M S E}^{2} = \frac{1}{N} \sum_{t = 1}^{N} {({G H I}_{t} - {\hat{G H I}}_{t})}^{2}$

${R M S E}^{2} = σ_{e}^{2} + e^{2}, e = \frac{1}{N} \sum_{t = 1}^{N} e_{t} (biais)$

$[{\hat{G H I}}_{t} - z_{1 - α / 2} R M S E; {\hat{G H I}}_{t} + z_{1 - α / 2} R M S E]$

$G H I$ (with a bar) is the mean of the measurements, $\hat{G H I}$ (with a hat) the forecast, N their number, $e_{t} = {G H I}_{t} - {\hat{G H I}}_{t}$ the residual, $σ_{e}$ the residual standard deviation, ē their mean (the bias), $z_{1 - α / 2}$ the quantile of the standard normal distribution (1.96 for 95%, see the normal distribution). This band has the same width at every hour: that is the homoscedasticity assumption. It contains 1 − α of the measurements only if the residuals follow a normal distribution.

Conformal prediction, without leakage

Split conformal prediction assumes no distribution. It measures directly, on data the model has never seen, how wrong it is, and derives the margin from that.

Training, 1998-2005: the model coefficients are estimated.
Calibration, 2006: this year is forecast without touching the model, and the 4380 daytime absolute errors |e| are sorted. The margin q is the one that leaves at most 5% of the errors above it: here q = 174 W/m² for AR at 95%.
Test, 2007: the interval is $\hat{G H I} \pm q$ , and the measurements inside are counted. Nothing from 2007 was used before.

Why three separate years? Calibrating on the training years would use errors the model has already “seen”, smaller than the ones it will really make: the margin would be too narrow. That is information leakage. Likewise, the AR order and the ELM size were chosen on 2005 only.

Result for AR at 95%: 95.2% coverage in 2007, with a mean band width of 318 W/m².

Three ways to set the margin

The three versions in the “Prediction interval” menu follow the same idea, that of our papers on the role of the clear-sky model (Solar Energy, 2025, Appendix C) and on complex-valued forecasting (JRSE, 2022). The interval is written $\hat{G H I} \pm k_{α} \hat{σ}$ . The multiplier $k_{α}$ is not taken from the normal distribution: it is read from a lookup table built from the errors actually made. What changes from one version to the next is the number of tables: one, one per hour, or one per clear-sky level. On this page, each table is computed on the year 2006, apart from training, to avoid leakage.

AR, nominal level 95%, test 2007. Coverage (PICP), mean band width and CRPS (score of the whole distribution, lower is better, see below).
Margin	Number of tables?	PICP 2007?	Clear days?	Variable days	Overcast days	By hour, 8:00–17:00?	Band width?	CRPS?
Gauss ±1.96 RMSE (reference)	0 (normal law)	93.4%	98.4%	85.2%	89.7%	86–98%	280 W/m²	39.0 W/m²
1. A single margin	1	95.2%	98.9%	88.5%	94.0%	90–98%	318 W/m²	37.4 W/m²
2. One margin per hour	14	94.5%	97.6%	89.2%	92.7%	93–98%	319 W/m²	37.3 W/m²
3. One margin per clear-sky class?	7	94.6%	98.1%	89.7%	90.2%	92–97%	325 W/m²	37.2 W/m²

1. A single margin: the lookup table

This is the method of Appendix C of the 2025 paper. The forecast errors are computed, then, for each level 1 − α, the multiplier $k_{α}$ giving the right coverage is found. It is stored in a table, which is then used for every new forecast. In the paper, the table depends on the site and the horizon; here, one site and one horizon, so a single table. With $\hat{σ} = R M S E$ , the margin is $q = k_{α} R M S E$ , the quantile of the absolute errors: this is exactly split conformal prediction.

The table already shows why the normal distribution goes wrong. On the 2006 errors:

Level 1 − α	80%	90%	95%	99%
$k_{α}$ read from the errors?	0.96	1.60	2.25	3.77
$z_{1 - α / 2}$ of the normal law?	1.28	1.64	1.96	2.58

At 80%, the errors call for a smaller multiplier than the normal law: many hours have an error close to zero (the peak of the histogram). At 99%, a much larger one is needed: large errors, when a cloud passes, are more frequent than a Gaussian bell predicts (the tails). A single table corrects the overall level; it keeps the same margin at every hour.

2. One margin per hour

This variant is not in our papers: it is an adaptation proposed for this page. Instead of one table, one is built per hour of the day (14 hours have Sun often enough in 2006). At 9:00, the forecast is only compared with the 9:00 errors of the calibration year. The margin then follows the Sun: from 69 W/m² at the extreme hours to 222 W/m² around noon, for AR at 95%. The hour is known before the measurement, so it may be used.

Result: from 8:00 to 17:00, coverage ranges from 93 to 98%, instead of 86 to 98% for the Gaussian interval, for an almost identical mean band width. At sunrise and sunset hours, which only have Sun part of the year, the table rests on few errors and coverage stays unstable (from 84 to 99%). In statistics, running a separate conformal computation in each group is called “Mondrian” conformal prediction.

3. One margin per clear-sky class

The 2022 paper makes the interval proportional to the clear-sky irradiance: $\hat{G H I} \pm μ {G H I}_{c s} \hat{σ}$ . The idea is the same here, in a more robust form. ${G H I}_{c s}$ is the irradiance that would be measured at that hour without any cloud (Solis model): it can be computed in advance for any date, from the Sun’s position and an average atmosphere only. The hours are sorted into 7 classes of ${G H I}_{c s}$ , slices of 150 W/m² (0–150, 150–300 … 750–900, above 900 W/m²), and one table is built per class. Its advantage over the hour: the class accounts for the season (noon in December resembles 9:00 in June). The margin ranges from 72 to 266 W/m².

This is not a cloudiness class. The class says how much Sun there could be, not whether there will be clouds, which is not known in advance. And slicing by 150 W/m² is an arbitrary choice, made before looking at 2007, which leaves at least 168 calibration errors in each class. This choice does not drive the result: with slices of 100, 200 or 300 W/m², or with 4, 7 or 10 equal-count classes, the 2007 PICP of AR at 95% stays between 94.6 and 94.8%, the band width between 321 and 328 W/m², and the CRPS between 37.2 and 37.4 W/m².

Why classes rather than a proportional margin? Dividing the error by ${G H I}_{c s}$ gives huge ratios at low-Sun hours, where ${G H I}_{c s}$ is tiny, and the interval widens everywhere: 682 W/m² mean band width for AR at 95% (94.5% coverage), against 325 W/m² with classes.

What none of the three fixes

Variable days stay below target (88.5, 89.2 and 89.7%), and clear days above. Neither the hour nor the clear sky tells whether a cloud is coming. That requires a quantity known in advance that measures how agitated the sky is: the role of the volatility $\hat{σ}$ forecast by complex-valued forecasting. The next step would therefore be the full interval of the 2022 paper, $\hat{G H I} \pm k_{α} {G H I}_{c s} \hat{σ}$ , with its table calibrated on a separate year.

Formulas: conformal quantile, lookup table and per-group margins

$s_{i} = | {G H I}_{i} - {\hat{G H I}}_{i} |, i = 1, \dots, n$

$q = s_{(k)}, k = ⌈ (n + 1) (1 - α) ⌉$

$[max ({\hat{G H I}}_{t} - q, 0); {\hat{G H I}}_{t} + q]$

$s_{(k)}$ is the k-th smallest of the n values $s_{i}$ , and ⌈·⌉ rounds up to the next integer; the “+ 1” corrects for the finite sample size. If calibration and test errors are exchangeable (same distribution, order irrelevant), the probability that a new measurement falls in the interval is at least 1 − α, over all hours. It is not a guarantee hour by hour. Hourly data are neither independent nor exchangeable (one hour resembles the previous one, and 2007 is not 2006): this theoretical guarantee therefore no longer applies automatically. Coverage is measured on 2007, without inferring a guarantee from it: 95.2% over all hours, 88.5% on variable days. The first line below links q to the multiplier $k_{α}$ of the lookup table. Versions by hour and by clear-sky class: the same computation is made inside each group g, with that group’s errors only, and the interval at hour t uses the margin $q_{g (t)}$ of its group, where g(t) is the hour of t or the class of ${G H I}_{c s, t}$ .

$q = k_{α} R M S E ⟺ k_{α} = \frac{q}{R M S E}$

$[max ({\hat{G H I}}_{t} - q_{g (t)}, 0); {\hat{G H I}}_{t} + q_{g (t)}]$

The lower bound is set to 0 because irradiance is never negative; this does not change the coverage.

PICP and mean band width

PICP: Prediction interval coverage probability: the share of measurements that fell inside their interval, in %. It is compared with the nominal level: at 95%, a PICP of 95% is reliable, 90% too optimistic, 99% too cautious.
Mean band width: The band width is the distance between the upper and lower bounds of the prediction interval, in W/m². Its mean over daytime hours is given; in brackets, the same relative to the mean daytime GHI (372 W/m²). Called MIL (mean interval length) in our papers.

Read both together. A huge interval always meets its target and is useless; a narrow interval that misses it is misleading. The good interval meets the target with the smallest width. These scores only count daytime hours: at night, forecast and measurement are 0.

Formulas: PICP and band width

$P I C P = \frac{100}{N} \sum_{t = 1}^{N} 1 [L_{t} \leq {G H I}_{t} \leq U_{t}]$

$W = \frac{1}{N} \sum_{t = 1}^{N} (U_{t} - L_{t})$

$L_{t}$ and $U_{t}$ are the lower and upper bounds, N the number of daytime hours, 1[·] equals 1 if the condition is true and 0 otherwise, $W$ is the mean band width.

The CRPS: scoring the whole distribution at once

PICP and band width judge an interval at a single level. Yet each method in the menu actually provides a whole family of nested intervals, from 1% to 99%: a predictive distribution, the probability law expected for the next measurement. The CRPS (continuous ranked probability score) scores this distribution in one go. It is expressed in W/m², like an error, and lower is better. It penalises both a badly placed distribution (the measurement falls where it was judged unlikely) and an overly spread one.

For a forecast without a margin, the CRPS is exactly the absolute error: its mean is the MAE (mean absolute error). This is why our 2022 paper presents it as the probabilistic counterpart of the MAE. The comparison is therefore direct: the gap between MAE and CRPS measures what the margin brings. The nCRPS divides the CRPS by the mean daytime GHI (372 W/m²), to compare sites or periods; the tool shows it next to the CRPS.

Daytime CRPS in 2007 (W/m²), by model and way of setting the margin. MAE: the forecast alone, without a margin.
Model	MAE?	CRPS, Gauss	A single margin	By hour	By clear-sky class
Persistence	110.6	75.2	74.6	71.4	73.7
Smart persistence	47.0	40.0	38.2	38.1	38.0
AR	46.7	39.0	37.4	37.3	37.2
ELM	44.9	37.5	35.9	35.8	35.8

For AR, the margin takes the error from 46.7 W/m² (MAE) to a CRPS of 37.2 W/m², 20% less. For every model, margins read from the errors do better than the normal law. For AR, smart persistence and ELM, the three conformal versions are almost tied: conditioning on the hour or the clear sky spreads the coverage better but changes the overall score little. Persistence is the exception: the per-hour margin clearly improves it (74.6 → 71.4 W/m²), because its errors depend strongly on the hour (it ignores that the Sun rises or sets). Finally, the models rank as for the MAE: the CRPS depends first on the quality of the central forecast.

The CRPS does not depend on the nominal level chosen in the tool, since it uses all levels at once. The computation follows Appendix D of our 2025 paper, which describes the distribution by its quantiles rather than by a formula. Here, 99 percentiles: percentile τ (the value the measurement has probability τ of not exceeding) is a bound of the interval of level $| 2 τ - 1 |$ . For example, percentiles 2.5% and 97.5% are the bounds of the 95% interval.

Formulas: CRPS, percentiles and nCRPS

${C R P S}_{t} = \int_{- \infty}^{+ \infty} {(F_{t} (x) - 1 [x \geq {G H I}_{t}])}^{2} d x$

${C R P S}_{t} \approx \frac{2}{M} \sum_{i = 1}^{M} ρ_{τ_{i}} ({G H I}_{t} - Q_{τ_{i}, t}), τ_{i} = \frac{i}{100}$

$ρ_{τ} (u) = u (τ - 1 [u < 0])$

$C R P S = \frac{1}{N} \sum_{t = 1}^{N} {C R P S}_{t}, n C R P S = \frac{C R P S}{G H I}$

$F_{t}$ is the predicted cumulative distribution function for hour t (the predicted probability that the measurement is below x), $Q_{τ, t}$ its percentile τ, M = 99 the number of percentiles, and $ρ_{τ}$ the “pinball loss”, which weighs more a measurement falling far into a tail judged unlikely. The second line approximates the first from the percentiles; the gap vanishes as their number grows. If all percentiles equal the forecast, the second line gives $| {G H I}_{t} - {\hat{G H I}}_{t} |$ : on average, the MAE. $G H I$ is the mean GHI of the N daytime hours.

The four models

All forecast the next hour (1 h horizon) from past measurements. The rules are those of the course: zero forecast when the Sun is below the horizon, never a negative forecast.

Persistence: The next hour will be like the last one: ${\hat{G H I}}_{t + 1} = {G H I}_{t}$ . Simple, but it ignores that the Sun rises or sets.
Smart persistence: The gap to the clear sky is kept: ${\hat{G H I}}_{t + 1} = {G H I}_{t} - {G H I}_{c s, t} + {G H I}_{c s, t + 1}$ , where ${G H I}_{c s}$ is the clear-sky (cloudless) irradiance computed by the Solis model.
AR: Autoregression: a weighted sum of the last 2 values of the additive clear-sky index (GHI − ${G H I}_{c s}$ , centred), coefficients estimated by least squares.
ELM: Extreme learning machine: a small neural network (20 hidden neurons drawn at random then fixed, only the output weights are computed), on the same inputs as AR. It is explained on the frugal forecasting page. Its CI is computed with the hidden neurons held fixed.

Daytime errors on 2006 (RMSE): persistence 132 W/m², smart persistence 80 W/m², AR 78 W/m², ELM 76 W/m². Comparing these models against horizon: the forecast comparison tool.

Formulas: models and confidence interval

$y_{t} = {G H I}_{t} - {G H I}_{c s, t} - μ$

${\hat{y}}_{t + 1} = \sum_{j = 1}^{p} w_{j} y_{t - j + 1}$

${\hat{G H I}}_{t + 1} = max ({\hat{y}}_{t + 1} + μ + {G H I}_{c s, t + 1}, 0)$

${\hat{G H I}}_{t + 1} \pm z_{1 - α / 2} s \sqrt{𝐱_{t}^{⊤} {(𝐗^{⊤} 𝐗)}^{- 1} 𝐱_{t}}$

${(𝐗^{⊤} 𝐗)}^{- 1} \to {(𝐇^{⊤} 𝐇 + λ 𝐈)}^{- 1} 𝐇^{⊤} 𝐇 {(𝐇^{⊤} 𝐇 + λ 𝐈)}^{- 1}$

$\hat{y}$ is the forecast of y, μ the mean of GHI − ${G H I}_{c s}$ over training, w the coefficients, p = 2 the order, $x_{t}$ the vector of the last p values of y, X the matrix of all these vectors over training, s the standard deviation of the training residuals. For ELM (second line), x is replaced by h, the hidden-layer output, and (X^TX)⁻¹ by its regularized form, where H gathers the hidden outputs over training, I is the identity matrix and λ = 10⁻³. These formulas assume independent residuals: the real CI is somewhat wider, but stays of the order of 1 W/m².

Data and split

Hourly GHI measured at Ajaccio (Corsica), course file (station 20004002), converted from J/cm² to Wh/m²; each value is the mean of the past hour, labelled by its end hour. Years of 8,760 hours. Training 1998-2005, calibration 2006, test 2007. Types of day in 2007, from $k_{d}$ , the ratio of measured to clear-sky irradiance summed over the day: clear if $k_{d}$ ≥ 0.9 (205 days), variable if 0.5 ≤ $k_{d}$ < 0.9 (111 days), overcast if $k_{d}$ < 0.5 (49 days). The type of day is only known afterwards: it is used to judge, not to build the interval.

Keep in mind: one site, one horizon (1 h), hourly data; the nominal level can be set from 50 to 99%. The theoretical conformal guarantee assumes exchangeable data, which these hourly measurements are not, and would in any case only hold over all hours: the coverage shown is measured on 2007, not guaranteed. Results would change for another climate or another horizon.

Going further: complex-valued forecasting: set α and the interval width · the probabilistic animation and its lookup table · judge a forecast with NICE · variability and forecastability, sCV and F · courses

References: where the method comes from

The intervals on this page apply the method of these two papers:

Cyril Voyant, Milan Despotovic, Gilles Notton, Yves-Marie Saint-Drenan, Mohammed Asloune, Luis Garcia-Gutierrez. On the importance of clearsky model in short-term solar radiation forecasting. Solar Energy 294, 113490 (2025). doi:10.1016/j.solener.2025.113490 · arXiv:2503.07647
Interval $\hat{y} \pm k_{α} \hat{σ}$ with a multiplier $k_{α}$ read from a lookup table built on the observed errors, with no assumed distribution (Appendix C); PICP, MIL, interval score and CRPS computed from quantiles (Section 2.3 and Appendix D, used in the “CRPS” section). This is the “single margin” method of this page.
Cyril Voyant, Philippe Lauret, Gilles Notton, Jean-Laurent Duchaud, Luis Garcia-Gutierrez, Ghjuvan Antone Faggianelli. Complex-valued time series based solar irradiance forecast. Journal of Renewable and Sustainable Energy 14(6), 066502 (2022). doi:10.1063/5.0128131 · arXiv:2203.06701
Interval proportional to the clear sky, $\hat{G H I} \pm μ {G H I}_{c s} \hat{σ}$ , with a multiplier μ calibrated on data rather than taken from the normal law, and a forecast volatility $\hat{σ}$ . Basis of the “by clear-sky class” version and of the proposed next step. The CRPS is presented there as the probabilistic counterpart of the MAE (Figure 4).

What is specific to this page: calibration on a separate year (2006) rather than on training data, and the versions by hour and by clear-sky class. The general statistical framework, conformal prediction, and its per-group version, called Mondrian, are described by V. Vovk, A. Gammerman and G. Shafer, Algorithmic Learning in a Random World, Springer, 2005 (doi:10.1007/b106715).

Cite this page

APA

Voyant, C. (2026). Probabilistic solar forecasting: PICP and conformal intervals. https://www.cyrilvoyant.com/tools/probabilistic-solar-forecasting/

BibTeX — For a LaTeX report or thesis: the page, then the two papers whose method it follows.

@misc{voyant2026probabilistic,
  author = {Voyant, Cyril},
  title = {Probabilistic solar forecasting: PICP and conformal intervals},
  year = {2026},
  howpublished = {\url{https://www.cyrilvoyant.com/tools/probabilistic-solar-forecasting/}},
  note = {Accessed [date]}
}

@article{Voyant2025Importance,
  author = {Voyant, Cyril and Despotovic, Milan and Notton, Gilles and Saint-Drenan, Yves-Marie and Asloune, Mohammed and Garcia-Gutierrez, Luis},
  title = {On the importance of clearsky model in short-term solar radiation forecasting},
  journal = {Solar Energy},
  year = {2025},
  volume = {294},
  pages = {113490},
  doi = {10.1016/j.solener.2025.113490}
}

@article{Voyant2022Complex,
  author = {Voyant, Cyril and Lauret, Philippe and Notton, Gilles and Duchaud, Jean-Laurent and Garcia-Gutierrez, Luis and Faggianelli, Ghjuvan Antone},
  title = {Complex-valued time series based solar irradiance forecast},
  journal = {Journal of Renewable and Sustainable Energy},
  year = {2022},
  volume = {14},
  number = {6},
  pages = {066502},
  doi = {10.1063/5.0128131}
}

Course archive (all versions) · Solar forecasting course (10.5281/zenodo.23620101).