Mathieu David

Faly H Ramahatana Andriamasomanana

P J Trombe

Philippe Lauret

F Ramahatana

Probabilistic forecasting of the solar irradiance with recursive ARMA and GARCH models

Keywords: probabilistic solar forecasts, clear sky index, ARMA, GARCH, operational framework, recursive least square

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Forecasts of the power output of solar renewables are required to improve their penetration rate into electricity grids and also to ensure the security of the supply-demand balance. Indeed, accurate forecasts allow a better scheduling of the energy resources and a better operation of the units commitment. The recent development of grid-connected storages associated with intermittent renewables (solar, wind and wave) also requires power forecasts in order to optimize their operational management [START_REF] Hernández-Torres | Technico-economical analysis of a hybrid wave power-air compression storage system[END_REF] [START_REF] Haessig | Energy storage sizing for wind power: impact of the autocorrelation of day-ahead forecast errors: Energy storage sizing for wind power[END_REF] [START_REF] Hanna | Energy dispatch schedule optimization for demand charge reduction using a photovoltaic-battery storage system with solar forecasting[END_REF]. In the case of solar renewables, the power output is directly related to the level of received solar irradiance. Thus, forecasting the solar irradiance is the key feature of the PV power prediction.

Many tools and methods for forecasting solar irradiance were developed for various time horizons and spatial resolutions [START_REF] Diagne | Review of solar irradiance forecasting methods and a proposition for small-scale insular grids[END_REF]. For day ahead forecasting, the most accurate methods are based on Numerical Weather Predictions (NWP) [START_REF] Perez | Comparison of numerical weather prediction solar irradiance forecasts in the US, Canada and Europe[END_REF]. For horizons of several hours, satellite based models are the most suitable [START_REF] &type=pdf Kühnert | Satellite-Based Irradiance and Power Forecasting for the German Energy Market[END_REF] [START_REF] Dambreville | Very short term forecasting of the Global Horizontal Irradiance using a spatio-temporal autoregressive model[END_REF]. And finally, for very short term horizons, from few minutes up to 2 hours, the literature is dominated by statistical approaches based on time series modeling [START_REF] Reikard | Predicting solar radiation at high resolutions: A comparison of time series forecasts[END_REF] [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF]. Most of these works are dedicated to the development of point forecasts that is the mean solar irradiance over the considered time step. No forecast is perfect. Prediction intervals are necessary to assess the uncertainty associated to the point forecasts. In the field of wind power prediction, an extensive literature exists about probabilistic forecasting [START_REF] Zhang | Review on probabilistic forecasting of wind power generation[END_REF] [START_REF] Jung | Current status and future advances for wind speed and power forecasting[END_REF]. But, in the solar energy domain, only few works were published recently in relation to probabilistic forecasts. A first approach used numerical weather predictions (NWPs) models in order to derive the uncertainty of the solar forecasts. These models provide forecasts for a large range of time horizons, from 1 hour up to 72 hours. [START_REF] Iversen | Probabilistic forecasts of solar irradiance using stochastic differential equations: PROBABILISTIC SOLAR FORECASTS BY SDES[END_REF] proposed a stochastic differential equation framework with NWP for modeling the uncertainty associated with the point irradiance forecasts. [START_REF] Alessandrini | An analog ensemble for short-term probabilistic solar power forecast[END_REF] assessed uncertainty of power forecasts of solar plants with an analog ensemble method. The uncertainty of the forecasts is obtained by gathering past conditions corresponding to analog sets of a reduced number of parameters provided by NWPs. [START_REF] Mathiesen | Geostrophic Wind Dependent Probabilistic Irradiance Forecasts for Coastal California[END_REF] proposed a postprocessing of different NWPs to obtain an ensemble probabilistic forecasting of the solar irradiance. The forecast system is based on studying the correlation of uncertainty to local meteorological conditions describing synoptic-scale atmospheric flow. Another simple approach has been proposed by [START_REF] Lorenz | Irradiance Forecasting for the Power Prediction of Grid-Connected Photovoltaic Systems[END_REF]. A normal distribution with zero mean and a standard deviation dependent on the solar zenith angle and the cloud situation is assumed to describe forecast errors of solar irradiance. A last approach is based on the probabilistic behavior of time series. [START_REF] Bacher | Online short-term solar power forecasting[END_REF] used a weighted quantile regression conditioned to a normalized power to produce the uncertainty associated to the point forecasts obtained with an autoregressive model (ARX). They applied this model to forecast the power output of a small solar plant for horizon ranging from 1 hour to 1 day. [START_REF] Bracale | A Bayesian Method for Short-Term Probabilistic Forecasting of Photovoltaic Generation in Smart Grid Operation and Control[END_REF] propose 1 hour ahead probabilistic forecasts of the clearness index. The method uses only ground measurements and it is based on a Bayesian autoregressive model associated with a Monte Carlo simulation to generate the density probability function of the error. [START_REF] Bessa | Probabilistic solar power forecasting in smart grids using distributed information[END_REF] proposed an original method to provide 6 hours ahead probabilistic forecasts of the output power of a small fleet of PV systems. This last method is based on a vector autoregressive framework (VAR). It enables to assess the uncertainty of the forecast from a spatially distributed network of measurements. This paper focuses on the very short term solar forecasting, from 10 minutes to 6 hours, applied to a single location. For these short horizons, the operational models should have low computational requirements in order to operate fast and the inputs must be available in real-time. The same issues are usually met in the financial domain. Furthermore, time series of solar irradiance recorded at high frequency (1-10 minutes) are composed of periods of low variability (clear sky or overcast sky) and periods with high variability (passing clouds). This behavior is also observed in financial time series that exhibit stable periods and clusters of volatility. This characteristic of times series is called heteroskedasticity [START_REF] Engle | Autoregressive Conditional Heteroscedasticity with Estimates of the Variance of United Kingdom Inflation[END_REF].

The first aim of this work is to test models commonly used in econometrics. The point forecasts are generated with the well-known AutoRegressive Moving Average model (ARMA) [START_REF] Tsay | Analysis of financial time series[END_REF]. The associated prediction interval is produced with a Genearalized AutoRegressive Conditional Heteroskedasticity model (GARCH) [START_REF] Bollerslev | Generalized autoregressive conditional heteroskedasticity[END_REF]. The combination of these two models has been extensively tested in econometrics but to our best knowledge it is the first time in the field of solar forecasting. The second objective is to propose an operational framework of the forecasting methodology. More precisely only past data of the measured solar irradiance are used to run the models and the parameters of the two models are estimated with a recursive method [START_REF] Ljung | Theory and Practice of Recursive Estimation[END_REF].

The remainder of this paper is organized as follows. Section 2 describes the data used to build and to test the different models. Section 3 discusses the errors metrics used to assess the accuracy of the point and probabilistic forecasts. Section 4 depicts the structure of the different forecasting models with a special emphasis on the recursive estimation of the parameters. Section 5 details the results of point forecasts while section 6 gives the results of the probabilistic forecasts. Finally, Section 7 gives some concluding remarks. In this work, we used data recorded at 6 different locations in the world (see Table 1). These sites were selected because they offered high frequency measurements of global horizontal solar irradiance (GHI) for long periods. They also represent a large variety of climates. The time step of record is different between the considered stations but never exceed 1 minute. For this study we computed 10-min. and 1hour averages of GHI directly from the raw data for 2 consecutive years. To get continuous and workable time series of GHI, we applied the following rules to treat the data:

Data

-To avoid the side effects induced by the low accuracy of the solar measurements at sunrise and sunset, we remove all the data for a solar zenith angle superior to 80°. Thus, the obtained time series do not contain the null night values. This filtering method removes less than 1% of the annual total solar energy. -If a gap smaller than 20 minutes is detected, a linear interpolation is done in order to fill the missing data with the initial time step. -If a gap longer than 20 minutes is detected, the whole day is removed from the time series.

The resulting data samples are continuous because the first data of the current day follow the last data of the previous day. Hence, the first forecasts of a day are done using the sky conditions observed before the sunset of the previous day. The measurements of the six selected stations are consistent and we detected only few missing data (table 1). The sites of Oahu and Tampon present a significant amount of missing days, around 15%, because their periods of record are shorter than two years.

A GHI time series can be divided into a deterministic part and a stochastic part. A clear sky model can accurately estimate the annual and diurnal cycles that represent the deterministic part [START_REF] Iqbal | An Introduction To Solar Radiation[END_REF]. The aim of this work is to forecast the stochastic part induced by the cloud cover. Therefore, in order to remove the deterministic part from the radiation time series, we compute the clear sky index (kt*) by dividing the measured GHI with the GHI observed under a clear sky (eq. 1). The resulting time series is almost stationary and a linear process, like an ARMA, can be used to model it.

(1)

!" * = !"# !"# !"#$%
In this work, we use the clear sky data provided by the McClear model [START_REF] Lefèvre | McClear: a new model estimating downwelling solar radiation at ground level in clear-sky conditions[END_REF] available for free on the SoDa website ("SoDa," 2015). This model use the AODs, water vapor and ozone date from the MACC project.

Figure 1 plots the clear sky index distribution and table 1 (last rows) gives the solar variability. These two indicators allow a better understanding of the characteristics of the solar irradiance for each site. In figure 1, one can notice a significant number of occurrences of overirradiance, that is to say a clear sky index above 1. These cloud enhancement events can be observed everywhere in the world [START_REF] Almeida | Extreme overirradiance events in São Paulo, Brazil[END_REF]. This phenomenon is more visible when using short time scales as in our study with 10-min mean of solar irradiance. The solar variability metric was defined by [START_REF] Hoff | Modeling PV fleet output variability[END_REF] as the standard deviation of the changes in the clear sky index (!∆!" ∆! * ). In our case we choose to compute the solar variability for a sampling intervals of 10 minutes and 1 hour that coincide with the time steps of our time series.

The continental stations of Desert Rock and Fort Peck experience a weather dominated by clear skies with a relatively low variability. For the 4 insular locations, the occurrences of clear skies are lower and the variability is higher than the continental sites. However some differences exist between them. For example, Oahu presents more occurrences of clear sky than Fouillole but these two locations have the same variability. Lauret et al. [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF] showed that the ability to generate an accurate point forecast with a linear process depends on the level of variability of the solar irradiance. Thus, the typology of the solar irradiance (i.e. radiation level, distribution and variability) of a site is a key factor to understand the level of error of a solar forecasting method.

Error metrics

Several attempts have been made to define the state-of-the art validation metrics for point forecasts [START_REF] Hoff | Modeling PV fleet output variability[END_REF]) [START_REF] David | Output and errors indicators for solar forecasting models[END_REF]) [START_REF] Coimbra | Overview of Solar-Forecasting Methods and a Metric for Accuracy Evaluation[END_REF]. In the realm of the solar forecasting community, the commonly used error metrics for point forecasts are: the Root Mean Square Error (RMSE), Mean Absolute Error (MAE), Mean Bias Error (MBE) and the Forecast Skill (FS). The formulas used to evaluate these metrics are given in appendix C. In their work, David et al. [START_REF] David | Output and errors indicators for solar forecasting models[END_REF] show that the RMSE and MAE metrics remain the best indicators. These absolute indices are usually normalized by the mean GHI of the test period as the utility industry wants to understand error in relative terms rather than in absolute terms [START_REF] Hoff | Modeling PV fleet output variability[END_REF]. In this work, in order to characterize the quality of the point forecasts, we provide the standard set of relative error metrics.

To evaluate and to rank the probabilistic forecasts, and following [START_REF] Iversen | Probabilistic forecasts of solar irradiance using stochastic differential equations: PROBABILISTIC SOLAR FORECASTS BY SDES[END_REF], (Alessandrini et al., 2015), [START_REF] Sperati | The "Weather Intelligence for Renewable Energies" Benchmarking Exercise on Short-Term Forecasting of Wind and Solar Power Generation[END_REF], [START_REF] Bessa | Probabilistic solar power forecasting in smart grids using distributed information[END_REF] and [START_REF] Pinson | Reliability diagrams for non-parametric density forecasts of continuous variables: Accounting for serial correlation[END_REF], we propose to use the reliability diagram, the rank histogram (or Talagrand histogram), the continuous ranked probability score (CRPS) and its associated skill score the CRPSS.

The reliability diagram is a well-known representation to verify a probabilistic forecast. We used the methodology defined by [START_REF] Pinson | Reliability diagrams for non-parametric density forecasts of continuous variables: Accounting for serial correlation[END_REF] to create it. The latter is specially designed for the case of continuous variables commonly met in the field of energy meteorology. On the horizontal axis, the reliability diagram gathers the forecasts into bins according to their cumulative probability. The frequency with which the meteorological parameter was observed to occur for each bin of the forecasts probability is then plotted against the vertical axis [START_REF] Hamill | Reliability Diagrams for Multicategory Probabilistic Forecasts[END_REF]. This representation is appealing because the deviations from perfect reliability (the diagonal) can be visually assessed [START_REF] Pinson | Reliability diagrams for non-parametric density forecasts of continuous variables: Accounting for serial correlation[END_REF].

The rank histogram is a tool for evaluating ensemble forecasts. They are useful for determining the reliability of probabilistic forecasts and for diagnosing errors in its mean and spread [START_REF] Hamill | Interpretation of Rank Histograms for Verifying Ensemble Forecasts[END_REF]. A rank histogram is generated by repeating the following two steps for every the time step of the verification sample: first sorting the members of the ensemble forecast and then tallying the rank of the corresponding observation. For a number of members N, the number of ranks of the histogram of an ensemble is N+1. A flat rank histogram shows that the members of an ensemble system are statistically indistinguishable from the observations [START_REF] Alessandrini | An analog ensemble for short-term probabilistic solar power forecast[END_REF]. A ∪-shape indicates a lack of spread while a ∩-shape indicates that the probabilistic forecasts are over-dispersive. In this work, we propose a parametric model to generate the probabilistic forecasts. To plot the rank histograms, we created ensemble forecasts based on this model. The members of the ensemble are generated with a random variable that follow the distribution law assumed in the parametric model.

The CRPS is also a commonly used verification index in order to compare a full probabilistic distribution with the observations, when both are represented as cumulative distribution functions (CDF) (Sperati et al., 2015). This index contains a sense of distance of how far the forecast was found from reality [START_REF] Hersbach | Decomposition of the Continuous Ranked Probability Score for Ensemble Prediction Systems[END_REF]. In addition, for a deterministic forecast, the CRPS is the mean absolute error (MAE) and has a clear interpretation. The formulation of the CRPS is:

(

2) CRPS = ! ! ! !"# ! ! -! !"# ! ! ! !! !! ! !!! !"
where ! !"# ! ! is the cumulative distribution function (CDF) of the forecast and ! !"# ! ! is the CDF of the observation for the i th ensemble forecast/observation pair. N is the number of pairs. The CRPS has the same dimension as the forecasted variable and a lower value (0 is the perfect score) of the CRPS corresponds to a better performance. The interested reader can refer to [START_REF] Hersbach | Decomposition of the Continuous Ranked Probability Score for Ensemble Prediction Systems[END_REF] that gives a complete and detailed presentation of how to asses the CRPS. As the metrics used to assess the accuracy of the deterministic forecasts, the relative counterpart of the CRPS is obtained by dividing the absolute value of this metric by the average GHI of the test period.

As proposed by [START_REF] Bessa | Probabilistic solar power forecasting in smart grids using distributed information[END_REF], the forecast skill of the proposed model is evaluated by computing the improvement over a reference model in terms of CRPS. This comparative metric will be called further in this work the Continuous Rank Probability Skill Score (CRPSS). The reference model used to derive the CRPSS is a persistence ensemble (see sub-section 4.1). The formulation of the CRPSS is given in appendix C.

Forecasting methods

In this section the forecasting models and the recursive method used to estimate their parameters are presented. To clarify the equations we used the following notation. A variable with a hat (^) is a forecast. h is the horizon of forecast also called the lead time. t denotes the time when the forecasts are generated. We choose to create one model per horizon of forecast in order to avoid the error propagations usually observed while running the same model iteratively. The models and the method of regression presented in the following sub-section are expressed taking into account this approach.

Reference models

Following [START_REF] Perez | Validation of short and medium term operational solar radiation forecasts in the US[END_REF], we propose to test our forecasting method against reference models like persistence, smart persistence and climatology.

The persistence model is expressed as follows:

(

3) !" * !!! = !" * !
This model assumes that the clear sky index for each time horizon 'h' only depends on the previous clear sky index values, which means that the sky conditions remain invariant between time 't' and time 't+h'.

The next model represents an easy way to improve persistence model and it is called Smart Persistence (Pers smart ). It consists in forecasting clear sky index for each time horizon 'h' as the mean of the 'h' previous clear sky values [START_REF] Perez | A new operational solar resource forecast model service for PV fleet simulation[END_REF]. Smart-persistence model is defined by the following equation:

(4) !" * !!! = !"#$ !" * ! , … , !" * !!!
We also propose the climatological mean model, which is independent of the forecast time horizon. More precisely, this model performs a constant forecast of the clear sky index that corresponds to its mean historical value:

(

5) !" * !!! = !"#$ !" *
During this work, we used the mean clear sky index of the first year of records in order to forecast the second year.

Finally, we propose to use a persistence ensemble as reference model [START_REF] Alessandrini | An analog ensemble for short-term probabilistic solar power forecast[END_REF] to make a comparison with our probabilistic forecasts. The members of the persistence ensemble are derived from the 10 last records of the clear sky index. Thus, the probability density function (PDF) of a forecast is equal the PDF obtained from the 10 most recent observations:

(6) !"# !" * !!! = !"# !" * ! , !" * !!! , ⋯ , !" * !!!

ARMA model

The AutoRegressive Moving Average model (ARMA) is one of the most used time series models in forecasting. In particular, it has been extensively studied in renewable energy forecasting and, owing to its parsimony, it has turned out to be a very tough competitor. Applications already include among others forecasting of solar irradiance [START_REF] Bacher | Online short-term solar power forecasting[END_REF]. A general formulation of an ARMA(p,q) model with p autoregressive (AR) terms and q moving average (MA) terms is given by Tsay et al. [START_REF] Tsay | Analysis of financial time series[END_REF]. Its application to the h-ahead forecast of the clear sky index is given by the following equation:

(7) !" * !!! = ! ! + ! ! • !" * !!!!! ! !!! + ! ! • ! !!!!! ! !!!
with h=1,2,… the forecast horizon. The error term ! is the difference between the previous forecasts and observations as defined in the following equation:

(8) ! !!! = !" * !!! -!" * !!!
More broadly, an ARMA model can be formulated in a canonical way as follows:

(9) !" * !!! = ! ! ! ! ! with ! ! = ! ! , ! ! , ⋯ , ! ! , ! ! , ⋯ ! ! the
vector of parameters to be estimated and

! ! ! = 1, !" * ! , ⋯ , !" * !!!!! , ! ! , ⋯ , ! !!!!! the vector of inputs.
There are several methods to estimate the vector of parameters θ, with the two most widely implemented being the least squares (LS) and the Maximum Likelihood Estimation (MLE) methods. Here, we chose a variation of the LS method, namely the Recursive Least Squares (RLS) method. This method offers the advantage of reducing the computational cost for estimating the model with its parameters being updated in real-time as new data become available. This contrasts with more intensive estimation methods operating on a sliding window with the estimation process being carried out at each time step. The RLS method is particularly useful in an operational context where forecast have to be timely delivered. In the RLS method, ! is time varying and, by convention, it is subsequently denoted ! ! . The RLS method is detailed in the following sub-section 4.4.

GARCH Model

The family of the ARCH (AutoRegressive Conditional Heteroskedasticity) models, introduced by Engle in the early 80's [START_REF] Engle | Autoregressive Conditional Heteroscedasticity with Estimates of the Variance of United Kingdom Inflation[END_REF], is used to model the volatility (i.e. the variance) of time series in the financial domain. These models are particularly suitable to describe changes in variance over the time. This family of models is commonly used in econometrics in order to forecast the realized volatility of high frequency data (Andersen, 2000)(McAleer and[START_REF] Mcaleer | Realized Volatility: A Review[END_REF] or to assess the volatility of the error of point forecast models like linear regressions (AR, ARMA, etc.) [START_REF] Bollerslev | Generalized autoregressive conditional heteroskedasticity[END_REF]. In this work, we applied this second approach to compute the prediction intervals associated with the point forecasts generated by a recursive ARMA model. [START_REF] Bollerslev | Generalized autoregressive conditional heteroskedasticity[END_REF] proposed a Generalized AutoRegressive Conditional Heteroskedasticity (GARCH) model, which provides a more parsimonious representation than the simple ARCH model in many applications. GARCH models express the conditional variance as a linear function of lagged squared error terms and also lagged conditional variance terms [START_REF] Taylor | Volatility forecasting with smooth transition exponential smoothing[END_REF]. The general formulation of a GARCH(p,q) model, with p error terms, q conditional variance terms and an horizon of forecast h, is given by the equations 10 and 11 [START_REF] Tsay | Analysis of financial time series[END_REF]:

(10) ! !!! = ! ! ! !!!
where ! is a random variable uniformly distributed with a zero mean and a unitary variance.

(11

) ! !!! ! = ! ! + ! ! ! !!!!! ! ! !!! + ! ! ! !!!!! ! ! !!! with the error term ! !!! = ! !!! -! !!! .
The basic idea is that the series is either serially uncorrelated or with minor lower order serial correlations, but it is a dependent series. GARCH models attempt to capture such dependence in the return series [START_REF] Tsay | Analysis of financial time series[END_REF]. Figure 2 shows a representation of the serial correlations of the time series of forecasting errors. These errors are obtained for the site of Saint-Pierre with a recursive ARMA(p,q) and a lead time of 10 minutes. Figure 2(a) plots the autocorrelation function (ACF) of the errors, which suggests no serial significant correlations. Figure 2(b) shows the ACF of the squared errors of forecast. This second plot suggests that the time series contains autocorrelation. Similar ACF plots and characteristics can be observed for other lead times and for the other sites. [START_REF] Bollerslev | Generalized autoregressive conditional heteroskedasticity[END_REF] shows that such time series can be modeled with a GARCH(1,1). Other models of heteroskedasticity exist, like the ARCH model. With a deeper analysis of the correlation in the data, they could also be set up with the same objectives [START_REF] Boland | Spatial-temporal forecasting of solar radiation[END_REF]. In order to define the confidence intervals of the point forecasts from the variance derived with the GARCH model, we assume a Gaussian distribution of the error. Thus, for each point forecast, the density of probability of the error corresponds to a normal law ! !" * !!! , ! !!! . This assumption is not rigorously supported by the data. Figure 3 shows the statistical distribution of the error of forecast for two sites that present significant differences in sky conditions. The distribution of the errors for Desert Rock presents a strong peak around the mean (kurtosis) and a lack of symmetry (skewness). It cannot be considered as gaussian. The distribution of the errors for Fouillole is closer to a normal law but there are also some important differences. Applied to the modeling of the forecasting error, equation 11 can be expressed with a canonical formulation as follows:

(12) ! !!! ! = ! ! ! ! ! with ! ! = ! ! , ! ! , ⋯ , ! ! , ! ! , ⋯ ! ! the
vector of parameters to be estimated and

! ! ! = 1, ! ! ! , ⋯ , ! !!!!! ! , ! ! ! , ⋯ , ! !!!!! !
the vector of inputs. The error term ! is the difference between the point forecast and the observation of the clear sky index at time t (eq. 13).

(13

) ! ! = !" * ! -!" * ! A special issue concerning the GARCH model is the estimation of the past values of the variance (! ! ! , …, ! !!!!! !
). The ARMA process generates only one point forecast per time step and it is not possible to assess the variance of past errors. Engle define ARCH processes as time series with non-constant variances conditional on the past, but constant unconditional variances [START_REF] Engle | Autoregressive Conditional Heteroscedasticity with Estimates of the Variance of United Kingdom Inflation[END_REF]. This unconditional variance ! ! is obtained by applying unconditional expectation to ! !!! ! in equations 11 or 12 [START_REF] González-Rivera | Smooth-Transition GARCH Models[END_REF]. The unconditional variance is a function of the parameters ! ! and ! ! of a GARCH model and it can be used to initialize the GARCH process (i.e. ! ! ! = ! ! in eq. 11) [START_REF] Nelson | Inequality Constraints in the Univariate GARCH Model[END_REF].

Equation 14 gives the unconditional variance for a GARCH(1,1) process:

! = ! ! !!! ! !! !
In this work, we build a GARCH(1,1) model as proposed by [START_REF] Bollerslev | Generalized autoregressive conditional heteroskedasticity[END_REF]. The estimation of the vector of parameters ! of the GARCH model is done with the same recursive method used for the ARMA model.

Recursive estimation of the parameters

We consider a time series ! ! that can be modeled by a linear regression process:

(15) ! ! = ! ! ! ! ! + ! !
where ! ! is the vector of parameters and ! ! the vector of inputs at time t. ! ! is a random variable with a normal law distribution. A forecast of the time series ! !!! can be done with the same linear regression model considering a horizon of forecast h corresponding to h time steps:

(16) ! !!! = ! ! ! ! ! ↔ ! ! = ! !!! ! ! !!!
We set up a recursive least square (RLS) to estimate the vector of parameters ! ! of this time series. For this work, the RLS methodology has been slightly modified to take into account the different time horizons h. For instance, the update step of the recursive process must be lagged in order to take into account that the forecasted ! ! was derived at time t-h. The adaptation of the RLS method integrates a different time lag between the current time and the forecasted time. In the equation 17 and 18, the gain matrix ! ! is assessed using a lagged input vector at time t-h. This change from the initial formulation takes into account that the forecast of the time series at time t was done at time t-h.

(

) ! ! = ! !!! + ! ! ! !!! ! ! -! !!! ! ! !!! (18) ! ! = ! ! ! !!! - ! !!! ! !!! ! !!! ! ! !!! !!! !!! ! ! !!! ! !!! 17
where ! ! is the gain matrix and ! is a forgetting factor. For this work, we use a forgetting factor ! = 0.999.

For the ARMA model, the RLS method was used in a classical way [START_REF] Ljung | Theory and Practice of Recursive Estimation[END_REF]) like [START_REF] Madsen | Time Series Analysis[END_REF] for the estimation of the parameters of their ARX model. The methodology directly minimizes the quadratic error between the observations ! ! and the forecasts ! ! . For the GARCH model, it is not possible to directly apply the RLS method. The prediction error is different from the GARCH residual. [START_REF] Kierkegaard | Estimating GARCH models using Recursive Methods[END_REF] transform the residual by making the transformation ! ! ! = ! ! ! in equation 17 in order to assess the GARCH parameters with a pseudo-linear recursive method. This transform works for any situations where the conditional variance is independent from the sign of the prediction error ! ! . We use this convenient transform to run our RLS method.

Even if the cyclical effects of the solar path are removed by the use of the clear sky index, the time series still present short-term trends (hurricanes, cold fronts, anticyclone, etc.). It seems important that the method used to assess the vector of parameters takes into account these short-term trends. The recursive methodology with forgetting is able to catch these short-term trends.

Machine learning techniques

Regarding the point forecasts, we also propose to compare the results obtained with the recursive ARMA with those obtained by nonlinear machine learning techniques (see section 5). In a previous work [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF], three machine learning techniques namely Artificial Neural Networks (NNs), Gaussian Processes (GPs) and Support Vector Machines (SVMs) were employed for intraday solar forecasting. These machine-learning techniques are supervised learning methods or data-driven approaches. As a consequence, these techniques rely on the information content embedded in the training data in order to produce forecasts on unseen data. More precisely, the models' parameters are determined with the help of n pairs of input and output examples contained in the training data. Once the model is fitted, the model can be evaluated on a test dataset. This second phase, called the generalization phase, consists of evaluating on the test dataset, the ability of the technique to generalize, that is to say, to give correct outputs when it is confronted with examples that were not seen during the training phase.

In our context, ! = x ! , y ! !!! ! represents the training dataset. The vector x ! contains the p past values of the clear sky index for training and y ! refer to the corresponding value of the clear sky index for the horizon h of interest. The column vector inputs for all n training cases can be aggregated in the so-called n×p design matrix X and the corresponding model's outputs (or targets) are collected in the vector y so we can write ! = X, y . Similarly, considering n * test cases, we have ! * = X * , y * for the test dataset.

These methods are briefly described below and the interested reader can refer to [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF] for details regarding the optimization of these types of models.

Neural networks (NNs)

A NN with d inputs, m hidden neurons and a single linear output unit defines a non-linear parameterized mapping from an input vector x to an output y. For our application, the prediction for an input test vector ! * is given by the following equation:

(19) !" * !!! = y x; w = w ! ! !!! f w !" ! !!! ! * + b ! + b !
Each of the m hidden units are related to the tangent hyperbolic function f x = e ! -e !! e ! + e !! . The parameter vector ! = w ! , w !" , b ! , b ! , which contains a set of weights w ! , w !" and two biases b ! , b ! , governs the non-linear mapping and is estimated during the training phase.

Gaussian Processes (GPs)

GPs are stated as a kernel-based method. Indeed, it can be shown [START_REF] Rasmussen | Gaussian processes for machine learning, Adaptive computation and machine learning[END_REF]) that, given n training samples, the prediction for an input test vector ! * can be seen in terms of a linear combination of n kernel functions; each one centered on a training point. Therefore, the forecasted clear sky index is given by Eq. ( 20):

(20) !" * !!! = α ! ! !!! ! ! ! ! , ! *
where ! ! denotes the squared exponential covariance function

! ! x ! , x ! = σ ! ! exp ! ! ! !! ! ! !! !
and ! ! is the ith input training vector. The coefficients α ! are estimated during the training phase. σ ! ! and l are called hyperparameters of the covariance function. These hyperparameters control the model complexity and can be learned (or optimized) from the training data at hand [START_REF] Rasmussen | Gaussian processes for machine learning, Adaptive computation and machine learning[END_REF].

Support vector machines (SVMs)

The support vector machine (SVM) is another kernel based machine learning technique used in classification tasks and regression problems [START_REF] Cortes | Support-Vector Networks[END_REF]. Support vector regression (SVR) is based on the application of support vector machines to regression problems [START_REF] Smola | A tutorial on support vector regression[END_REF]. This method has been successfully applied to time series forecasting tasks [START_REF] Muller | Predicting Time Series with Support Vector Machines[END_REF]. In a similar manner as for the GPs, the prediction calculated by a SVR machine for an input test case ! * is given by Eq. ( 21):

(21) !" * !!! = ! ! ! !!! ! !"# ! ! , ! * + ! k !"# denotes the radial basis covariance function k !"# x ! , x ! = exp -γ x ! -x ! with
hyperparameter γ and b is a bias parameter.

In the case of SVMs (unlike GPs), it must be noted that not all the training patterns participate to the preceding relationship. Indeed, a convenient setting of the optimization problem enables to obtain a sparse solution. The latter means that only some of the coefficients α ! will be nonzero. The examples that come with non-vanishing coefficients are called Support Vectors.

Results for the point forecasts

The 2 years of data of the six considered sites are split into two sets as follow: the first year is used as training sample and the second year is used as test sample. All the metrics provided in this section are derived from the test samples (i.e. the second years of data). In the case of the recursive estimation of the parameters of the ARMA models, a calibration or training period is not required. Thus, to select the best orders of the ARMA(p,q) models, we directly used the first years of data of every sites. We selected the pairs of AR and MA lags that presented the lowest RMSE for each site and each forecasting horizon. To find this minimum, we varied the lags in the range [1;10]. Table 2 gives the AR and MA orders resulting from the cross validation. For a 10-min. granularity and for the sites with a low variability (i.e. Desert Rock, Fort Peck and Saint-Pierre), the AR and MA orders present the same evolution. As shown by Table 2, the AR order increases and the MA order decreases when the forecasting horizon is increasing. Obviously, for these sites that present a relatively low variability in the 10-min time scale, very short-term forecasts (horizons inferior to 30 minutes) are mainly influenced by the trend of the change of the sky condition. For the sites with a higher variability and also for a granularity of 1 hour, no clear trend was seen in the values of the orders.

1 hour 1 1 1 1 1 1 1 1 1 1 1 1 2 hours 2 1 1 1 2 2 4 1 2 1 9 1 3 hours 1 1 1 1 3 2 1 2 7 5 8 2 4 hours 2 2 3 3 2 1 1 1 7 2 9 1 5 hours 1 1 1 1 5 1 5 2 6 1 8 1 6 hours 1 1 1 1 5 1 1 7 6 1 7 1
Table 2: Orders of the recursive ARMA(p,q) models for different lead times, granularities and sites

The detailed results of the simulations can be found in the tables A.1 to A.4 of the appendix A. These tables give the values of the error metrics for every site, horizon and granularity. The proposed method presents a very low bias (MBE), inferior to 2% for all sites. Figures 4 and5 show the relative RMSE of the point forecasts respectively for the 10-min. granularity and for the 1-hour granularity. In terms of quadratic errors, the accuracy of the forecasts decreases when the lead-time increases. The recursive ARMA models outperform the two reference models (i.e. smart persistence and climatology) for all the horizons and granularities. As mentioned above, a comparison is also done with the machine learning techniques briefly described in section 4.5. The Gaussian Process technique from this previous work was not used for the 10-min. granularity because it led to a very long computation time. In terms of quadratic errors, the recursive ARMA slightly outperforms the machine learning models (Gaussian Process, Support Vector Machine and Neural Network) for a granularity of 1 hour. But for a granularity of 10 minutes and for the other error metrics (see tables 5 and 6 in appendices), the accuracy of the recursive ARMA is practically the same as the machine learning methods. The main advantage of the proposed ARMA model compared to the machine learning methods comes from the simplicity provided by the recursive estimation of the parameters. Indeed, this method does not require training period and the parameters are estimated only with the last available records using a simple algorithm. Moreover, an ARMA model has only few parameters to estimate and it does not require complicated method to manage its complexity. Hence, it runs very fast and it is easier to set up for an operational use.

The accuracy of the forecasts also depends on the sky conditions observed on a site. There is a clear relationship between the indices of variability (!∆!" !"!"# * and !∆!" !"!"# *

) and the ability to forecast accurately the solar irradiance. The stations of Oahu, Fouillole and Tampon experience the most variable sky conditions and they exhibit the less accurate forecasts. On the contrary, Desert Rock has the lowest index of variability and the accuracy of the forecasts is the best for this site. The more variable the sky conditions are, the less predictable the solar irradiance is. This relationship has already been highlighted in recent works [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF] [START_REF] Voyant | Mean absolute log-return and solar radiation forecastability[END_REF]. [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF] considering all the sites. [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF] considering all the sites.

Results for the probabilistic forecasts

An uncertainty envelope is associated to the point forecasts thanks to the GARCH(1,1) model (Figure 6). This probabilistic information is derived from the forecasts of the standard deviation of the error between the point forecasts and the observations. To generate the uncertainty, we assume that the distribution of the error of forecast follows a normal law. To assess the ability of the model to reproduce the distribution of his own error, we first use the well-known reliability diagram. This representation allows to graphically assess the reliability for a set of probabilistic forecasts. An example of a reliability diagram is given in 7, the probabilistic forecasts do not follow exactly the ideal line. For example, 54% of the observed errors occurred when forecast probability was 60%. In this case, the model is over confident. The probability of the forecasts is higher than the probability of the observations. For all the horizons, the reliability of the ARMA-GARCH models is almost same. Around a cumulative probability of forecasts of 70%-80% the reliability of the model is changing. Below this threshold, the model is over confident. Inversely, above this threshold, the model is under confident. Strangely, the model based on stochastic differential equations developed by [START_REF] Iversen | Probabilistic forecasts of solar irradiance using stochastic differential equations: PROBABILISTIC SOLAR FORECASTS BY SDES[END_REF] shows a similar behavior. In our case, this deviation of the probabilistic forecasts to the ideal line is mainly due to the assumption done about the Gaussian distribution of the error. Figure 9 gives the effect of this assumption on the reliability of the probabilistic forecast. The distance between a Gaussian distribution and the distribution of the error of the forecasting method is measured with the help of the Kolmogorov-Smirnov index (KSI). The KSI for each site is given in the table at the right part of figure 9. The smaller the KSI is, the closer the distribution to a normal law is. The distribution of the forecasting error of the site of Desert Rock is the most distant from a normal law. The worst reliability is also observed for this site. Inversely, the site of Tampon presents the closest distribution to a normal law and the best reliability of our test samples. Thus, the assumption done about the statistical distribution of the error of the point forecasts affects directly the reliability of the probabilistic forecast done with our GARCH model. In our case, the sites that exhibit the lowest error for the point forecasts have the worst reliabilities. Thus, the higher the solar variability is, the better the reliability of the probabilistic forecasts is. To complete the analysis done above about the reliability, we plot rank histograms to assess the consistency of our method. Figures 10 and 11 show the rank histograms of our ARMA-GARCH models and of the persistence ensemble forecasts. The rank histogram is initially designed to assess the spread of an ensemble forecast system. In the case of a parametric method, like our approach, a rank histogram can be done by generating an ensemble of forecasts from a random variable that follow the selected distribution. As for the persistence ensemble, we choose to generate 10 members that follow our assumption of normal distribution of the error. The mean comes from the ARMA model and the variance comes from the GARCH model. With a U-shape, the rank histograms of the persistence ensemble show clearly a lack of spread in the forecasts. With a peak of population on the right side, the shape of the rank histograms of our ARMA-GARCH model shows an excess of variability with a negative bias. We can suppose that these bias results from the symmetry of the PDF of the normal law. Indeed, if the sky is clear at time t, at time t+h the probability that the sky will be "less clear" is higher than a probability to be "more clear". This reality is not supported by our assumption of a Gaussian distribution of the error of forecast. The last investigated error metrics are the CRPS and the associated skill score the CRPSS. The CRPScombines both the error of the point forecasts and the reliability of the probabilistic quantiles. Figure 12 (Tables B.1 and B.2 in appendix B) gives the CRPS of each site for all the granularities and all the horizons. As the CRPS is equivalent to the MAE for a deterministic forecast, its trend is strongly influenced by accuracy of the point forecasts. Again, we can observe that the accuracy of the forecasts is directly dependent to the solar variability of the site. Even if the reliability of the probabilistic forecasts differs between the different sites, it affects very weakly the values of the CRPS. In figure 13, the CRPSS gives the improvement of the ARMA-GARCH model compared to the persistence ensemble. Overall, the ARMA-GARCH model outperforms the persistence ensemble. This improvement is higher for the shortest horizons of forecast. There is also a difference of improvement between the different sites. Again, the variability of the observed sky conditions is the key factor of these differences. By definition, the persistence ensemble has a good accuracy when the sky conditions remain stable. Thus, the improvement of our model is better for the sites that experience a higher variability. The results obtained for the analysis of the probabilistic forecasts are coherent with the formulation and the behavior of a GARCH model. For instance, a GARCH process is efficient to model time series that tend to "return to the mean" after an event of strong volatility. This mean volatility is also the unconditional volatility of the time series. For our time series of clear sky indices, the mean volatility (or variance) is site dependent and it has a value around 0.2 -0.3. These mean values do not correspond to a stable sky condition (i.e. clear sky or overcast) that exhibits an almost null variance. Thus, for stable sky conditions, the GARCH process tends to overestimate the volatility of the error.

Conclusions and outlook

In order to provide very short-term forecasts of the solar irradiance, operational tools with low computational requirements must be developed. The approach proposed in this work is built on models used for the high frequency trading in the financial domain. The ARMA-GARCH methodology is an efficient combination of models used to generate very short-term point forecasts of solar irradiance with confidence intervals. Furthermore the recursive estimation of the parameters of these two models facilitates the implementation of these methods in an operational framework. In other words, this approach can be set up using only the real time monitoring of the solar power.

The recursive ARMA model offers a simple and operational tool to generate point forecasts. Using only the actual value of the solar irradiance, this technique outperforms other statistical models as applied in this work. The association of an ARMA model with a GARCH model allows assessing the prediction intervals associated with the point forecasts. This combination is an alternative to the ensemble approach commonly used for generating probabilistic forecasts. Furthermore, a recursive ARMA-GARCH needs lower computational requirement than an ensemble method based on NWP and it is more suitable for very short term forecasting.

In this work, a normal distribution of the error of the point forecasts has been assumed. In reality, the error of the point forecasts does not fit with a Gaussian law. Following this assumption, the model is over confident when the forecast probability is below 60%-70% and under confident above this threshold. An improvement of the method could be achieved by assuming a distribution law that fits better to the data (e.g. a Student law).

New tools are currently developed in order to provide real time monitoring of the sky conditions like fisheye cameras. The proposed framework can easily integrate exogenous inputs (e.g. through the use of an ARMAX model) in order to improve the quality of the point forecasts. Thus, the recursive ARMA-GARCH model is not only an operational model but it is also an upgradeable tool. (C.4) !" = 1 -!"#$ !"#$!"#$%&' !"#$% !"#$ !"#$%&'()" !"#! !" * * 100 (C.5) !"#$$ = 1 -!"#$ !"#$%&'()*+ !"#$% !"#$ !"#$%&'()" !"#!$%&! !"#! !" * *100

Figure 1 :

 1 Figure 1: Statistical distribution of the 10-min. clear sky index for the 6 considered locations

Figure 2 :

 2 Figure 2: Sample ACF of the 10-min. horizon errors of forecast derived from a recursive ARMA for the site of Saint-Pierre ; (a) ACF of the errors ; (b) ACF of the squared errors

Figure 3 :

 3 Figure 3: Statistical distribution of the errors of forecast (bar chart) derived from a recursive ARMA model with a horizon of 30 minutes for the sites of Fouillole (a), Desert Rock (b) and the best fit curve corresponding to a normal law (line).

Figure 4 :

 4 Figure 4: Relative RMSE for horizons of forecast from 10 minutes to 60 minutes and a granularity of 10 minutes. In (a), the dashed lines are the error of the recursive ARMA model for each individual site. The bold line is the average RMSE of the recursive ARMA model considering all the sites. (b) Comparison of the reference models, the recursive ARMA and the models proposed in[START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF] considering all the sites.

Figure 5 :

 5 Figure 5: Relative RMSE for horizons of forecast from 1 hour to 6 hours and a granularity of 1 hour. In (a) the dashed lines are the error of the recursive ARMA model for each individual site. The bold line is the average RMSE of the recursive ARMA model considering all the sites. (b) Comparison of the reference models, the recursive ARMA and the models proposed in[START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF] considering all the sites.

  figure 7 for probabilistic forecasts done for the site of Tampon. The dashed line represents a perfect assessment of the errors by the model itself. As shown in figure7, the probabilistic forecasts do not follow exactly the ideal line. For example, 54% of the observed errors occurred when forecast probability was 60%. In this case, the model is over confident. The probability of the forecasts is higher than the probability of the observations.

Figure 6 :

 6 Figure 6: Example of point forecasts with confidence intervals of the GHI for 2 days in Oahu generated with the recursive ARMA-GARCH model. The forecasting time horizon is 1 hour and the granularity is 10 minutes.

Figure 7 :Figure 8

 78 Figure 7: Reliability diagram of the probabilistic forecasts for Tampon with a forecasting time horizon of 1 hour and a granularity of 1 hour

Figure 8 :

 8 Figure 8: Reliability diagrams of the recursive ARMA-GARCH model considering all the sites. (a) Horizons of forecast from 10 minutes to 60 minutes and a granularity of 10 minutes. (b) Horizons of forecast from 1 hour to 6 hours and a granularity of 1 hour.

Figure 10 :

 10 Figure 10: Rank histograms computed over all the lead times of the ARMA-GARCH model (a) and of persistence ensemble (b) for a 10-min granularity.

Figure 11 :

 11 Figure 11: Rank histograms computed over all the lead times of the ARMA-GARCH model (a) and of persistence ensemble (b) for a 1-hour granularity.

Figure 12 :

 12 Figure 12: CRPS of the recursive ARMA-GARCH. The dashed lines are the CRPS for each individual site. The bold line is the average CRPS considering all the sites. (a) Horizons of forecast from 10 minutes to 60 minutes and a granularity of 10 minutes. (b) Horizons of forecast from 1 hour to 6 hours and a granularity of 1 hour.

Figure 13 :

 13 Figure 13: CRPSS of the recursive ARMA-GARCH. The dashed lines are the CRPSS for each individual site. The bold line is the average CRPSS considering all the sites. (a) Horizons of forecast from 10 minutes to 60 minutes and a granularity of 10 minutes. (b) Horizons of forecast from 1 hour to 6 hours and a granularity of 1 hour.

  #$%&'(.! -GHI !"#$%&"'.! ! !!!

Table 1 :

 1 Main characteristics of the solar measurements used in this work

  Reliability diagrams of the recursive ARMA-GARCH for a horizon of 30 minutes and a granularity of 10 minutes. The one-sample Kolmogorov-Smirnov index (KSI) for a correspondence with a normal law is given in the table for the error of forecasts on the kt*.

		KSI
	Desert Rock	0.37
	Fort Peck	0.33
	Fouillole	0.31
	Oahu	0.31
	Saint-Pierre	0.33
	Tampon	0.29
	Figure 9:
	Observed frequency ratio	t+1 hour t+2 hours t+3 hours
		t+4 hours
		t+5 hours
		t+6 hours
		Ideal
	Forecast probability

Table A .

 A 1: Relative error metrics of the point forecasts done for each site by the recursive ARMA models using a granularity of 10 minutes

				)		)		)		)		)		)
		Horizon	Desert Rock	(mean = 548W/m 2	Fort Peck	(mean = 391W/m 2	Fouillole	(mean = 510W/m 2	Ohau	(mean = 499W/m 2	Saint-Pierre	(mean = 538W/m 2	Tampon	(mean = 456W/m 2
		1 hour	-0.2%	-0.7%	-0.9%	1.5%	-0.4%	1.5%
		2 hours	-0.5%	-1.2%	-1.8%	2.1%	-1.0%	0.5%
	MBE	3 hours 4 hours	-1.0% -1.3%	-1.5% -1.5%	-2.1% -2.2%	2.2% 2.0%	-1.3% -1.8%	0.3% 0.1%
		5 hours	-1.4%	-1.6%	-2.1%	1.8%	-1.9%	-0.1%
		6 hours	-1.4%	-1.6%	-1.9%	2.0%	-1.9%	-0.4%
	RMSE

Table B .

 B 2: CRPS and CRPSS of the probabilistic forecasts done for each site by the recursive ARMA models using a granularity of 1 hour

Appendix C: Error metrics formula

Acknowledgement

The authors would like to thank the laboratory LARGE from the University of Guadeloupe, the National Renewable Energy Laboratory (NREL) and the meteorological network SURFRAD for providing their ground measurements. This work was supported by the French Environment and Energy Management Agency (ADEME) under the SOLFIN project.

-0.1% -0.2% -0.3% 0.5% -0.2% 0.4% 20 min.

-0.2% -0.3% -0.6% 0.6% -0.3% 0.7% 30 min.

-0.3% -0.5% -0.8% 1.0% -0.4% 1.0% 40 min.

-0.3% -0.6% -0.9% 1.3% -0.5% 1.1% 50 min.

-0.4% -0.7% -1.1% 1.3% -0.6% 1.2% 60 min.

-0.4% -0.8% -1.3% 1.4% -0.6% 1.2%

RMSE