Semiparametric prediction models for variables related with energy production

González-Manteiga et al. Journal of Mathematics in Industry (2018) 8:7 https://doi.org/10.1186/s13362-018-0049-0 RESEARCH Open Access Semiparametric prediction models for variables related with energy production Wenceslao González-Manteiga1,2†, Manuel Febrero-Bande1,2*†and María Piñeiro-Lamas3† *Correspondence: [email protected] 1MODESTYA group, Technological Institute for Industrial Mathematics (ITMATI), Santiago de Compostela, Spain 2Dept. of Statistics, Mathematical Analysis and Optimization, Fac. of Mathematics, Universidade de Santiago de Compostela, Santiago de Compostela, Spain Full list of author information is available at the end of the article †Equal contributors Abstract In this paper a review of semiparametric models developed throughout the years thanks to an extensive collaboration between the Department of Statistics and Operations Research of the University of Santiago de Compostela and a power station located in As Pontes (A Coruña, Spain) property of Endesa Generation, SA, is shown. In particular these models were used to predict the levels of sulphur dioxide in the environment of this power station with half an hour in advance. In this paper also a new multidimensional semiparametric model is considered. This model is a generalization of the previous models and takes into account the correlation structure of errors. Its behaviour is illustrated in a simulation study and with the prediction of the levels of two important pollution indicators in the environment of the power station: sulphur dioxide and nitrogen oxides. Keywords: Semiparametric prediction models; Pollution indicators; Cointegration 1 Introduction: an environmental problem The coal-fired power station in As Pontes is one of the production centers owned by Endesa Generation SA in the Iberian Peninsula. It is located in the town of As Pontes de García Rodríguez, northeast of A Coruña province. This power station was designed and built to make use of lignite from the mine located in its vicinity. This solid fuel was characterized by its high moisture and sulphurcontents anditslowcalorificvalue.Throughouttheyearstheplanthasundergoneseveraltransformation processes in their facilities with the aim of reducing emissions of sulphur dioxide (SO2). The power station completed its last adaptation in 2008 to consume, as primary fuel, imported subbituminous coal, characterized by its low sulphur and ash contents. Thelocationofthepowerplantclosetonaturalsitesofhighecologicalvalue,suchasthe NaturalPark AsFragasdoEumeandexistinglegislation,meanthatithasexistedsincethe beginning a great concern for its impact on the environment. Therefore the station has a Supplementary Control System of Air Quality that allows it to make changes in operating conditions in order to reduce emissions when the weather conditions are adverse to the spread of the emitted smoke plume,specifically containing SO2,andtherearesignificant episodesofimpairedairquality.Spanishlaw,byrulesandregulations,setsmaximumconcentrationsthatcanbeachievedforthesegasesinagivenperiodoftime.Inparticular,for this plant the only limit that might be exceeded at any time, is one that is established on ©The Author(s) 2018. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. González-Manteiga et al. Journal of Mathematics in Industry (2018) 8:7 Page 2 of 16 the hourly mean (continuously computed) from the concentration of SO2in the soil, the value of 350 μg/m3. The problem is to be able to predict, using the information received continuously at sampling stations and the past information, the future values for SO2levels. Statistical forecast models are the key to get these predictions and suggest a course of action to the plant operators. In recent years, new statistical models have been designed to obtain the simultaneous prediction of two pollution indicators in the environment due to the changes in the environmentallegislation,inthepowerstationitself,andtheconstructionofanewnaturalgas combined cycle station in the vicinity. The fuels that are going to be used make that the maininterestliesinpredictingthevaluesofthenitrogenoxides(NOx)whichisemittedby both facilities simultaneously with the values of SO2which is only emitted by the power station. All these changes have created a new problem: predicting hourly mean concentrations ofsulphurdioxideandnitrogenoxides,measuredintheenvironmentofthetwofacilities. Faced with this new approach, the statistical forecast models are again an effective tool. Thus, a multidimensional prediction general model is designed (see Sect. 3). 2 Methods: one-dimensional predictive models 2.1 Models designed to solve the environmental problem ResultingfromthecollaborationoverthepastyearsbetweentheDepartmentofStatistics and Operations Research at the University of Santiago de Compostela and the Environment Section of the power station, the Integrated System of Statistical Prediction of the Immision (SIPEI, in Spanish) have been created employing statistical models to provide predictions for the levels of SO2with a half an hour horizon. Due to data availability with minutal frequency in real-time and current legislation, the hourly mean is considered from both of the values of SO2and NOx, for predictions of future values of both pollutants. Thus, two time series are constructed, X1,tand X2,t,for whichthesubscripttrepresentsaminutalinstant,andeachvaluewillbeanaverageofthe actual values for the last hour: X1,t=1 60 59  i=0 SO2(t–i)andX2,t=1 60 59  i=0 NOx(t–i), where SO2(t)andNO x(t) represent the concentration of SO2and NOx,respectively,at time t,measuredinμg/m3. The series of hourly SO2means has a characteristic behaviour, highly influenced by weather conditions and local topography. It takes values close to zero for long periods of time, and it can suddenly and sharply increase (episodes) in bad meteorological conditions for the dispersion of the smoke plume. Nowadays, the series of hourly NOxmeans has a similar behaviour to that of SO2, but on a smaller scale (see Fig. 1). The main objectiveofthedevelopedstatisticalmodelsistopredicttheepisodes,soourinterestiscentred on the values that occur less frequently along the time series. Becauseofthis,akindofmemorycalledHistoricalMatrixwasdesigned(Prada-Sánchez and Febrero-Bande [14]), whichwillbeessentialtothebehaviour ofall developed models so far. This matrix is composed of a large number of vectors based on (Xt–l,...,Xt,Xt+k): González-Manteiga et al. Journal of Mathematics in Industry (2018) 8:7 Page 3 of 16 Figure 1 Episode depicted in one of sampling stations. The one hour mean of SO2and NOxare, respectively, drawn in red and orange realdataofbihourlySO2orNOxmeans,chosensoastocoverthefullrangeofvariablein question and make the role of historical memory. To ensure that cover the entire range of the variable, the matrix is divided into blocks according to the level of the response variable, Xt+k. To update the memory, in every instant, when a new observation is received, the historical matrix is renewed in the following way: the class to which the new observation belongs is found and then the oldest datum in such class leaves the matrix and the new observation enters it. With a sample built this way, makes sure that always have updated information on the full variation range of the interest variable, and over the years this concept has been adapted to the different statistical techniques used. 2.1.1 The first semiparametric model In the early years of development, the data transmission frequency to SIPEI was pentaminutal, and also, the legislation in force at that time established the limit values for the two hour mean of the SO2. For this reason, the prediction models for SO2levels initially worked with series of bihourly means. The objective was to obtain the prediction, with a half an hour horizon, for this time series. Therefore, each time it receives a new observation, Xt,ithastopredictthevalueatsixtimesahead,Xt+6. Asemiparametric approachwas considered (García-Jurado et al. [8]) which generalizes the traditional Box–Jenkins models as follows: Xt+κ=ϕκ(Xt,Xt–l)+Zt+κ,κ,l∈Z+, where Zthas an ARIMA structure of mean zero independent of Xt(Box et al. [1]). González-Manteiga et al. Journal of Mathematics in Industry (2018) 8:7 Page 4 of 16 In particular at each time t, the regression function ϕ6(Xt,Xt–1)=E(Xt+6/Xt,Xt–1)isestimatedwiththewell-known Nadaraya–Watson kerneltype estimator(seeNadaraya[13] and Watson [19]) using the information provided by the historical matrix. The second stepisto calculatetheresidualtimeseries ˆ Zt–64,...,ˆ Ztrelativetothe last sixhours, where ˆ Zi=Xi–ˆ E(Xi/Xi–6,Xi–7)foreachiandfitsanappropriateARIMAmodelforit.Finallywe get the Box–Jenkins prediction of ˆ Zt+6. The final point prediction proposed is given by: ˆ E(Xt+6/Xt,Xt–1)+ˆ Zt+6. 2.1.2 Partially linear model Theinformationused bytheprevious semiparametric modelstoobtainthepredictionsis thepastofthetimeseries;howeverit mightbeusefultointroduceadditionalinformation in order to improve these predictions. Specifically, meteorological and emission variables have been used with, the so-called partially linear models (Prada-Sánchez et al. [15]) to estimate bihourly mean values of SO2with one hour in advance. Dataintheformof(Vt,Zt,Yt)isconsidered,whereVtisavectorofexogenousvariables, Zt=(Xt,Xt–l)andYt=Xt+12 being Xtthe series of bihourly SO2means; and it is assumed thatthisseriesconformtothefollowingpartiallylinearmodel:Yt=Vt tβ+ϕ(Zt)+t,where tisanerrortermofmeanequalstozero. This model can easily estimated following Speckman [18] and allow us to extend the horizon to one hour maintaining the same level of accuracy as the semiparametric model for half an hour horizon. In any case, the incorporation of external information slightly improves theprediction because the measure point for the meteorological variables is located at 80 m over ground level which is relatively far away (and so, uncorrelated) respect to the typical height of the emitted smoke plume (above 800 m over ground level). Emissioninformationisalsooflittleinterestbecausethesesignalsarealmostconstantspecially when the facility is working not describing at all the reasons that make the smoke plume fallstotheground.Bythesereasons,meteorologicaloremissioninformationwasnotconsidered in the following models. 2.1.3 Neural networks The change in the interest series established by the European Council Directive 1999/30/CE, from bihourly means to hourly means, causes the time series to be less smooth. At the beginning, the previous semiparametric model was adapted to work on thenewseriesofhourlymeans.Theresultsshowedaconsiderableincreaseintermsofthe variability of the given predictions, regarding the results usually obtained for the series of twohourmeans. In an attempt to improve the response given by the SIPEI, and in particular, its point predictions with half an hour horizon, new predictors based on neural networks models were developed (Fernández de Castro et al. [6]). A neural network model has been designed to provide predictions of one hour mean values of SO2with half an hour in advance. It consists of an input layer, one hidden layer and an output layer. The number of nodes in the output layer is determined by the size of the response to be obtained from the network; in this case interested in a prediction for Xt+6.Asinputtothenetworkithasbeentakenthebidimensionalvector (Xt–3,Xt)andthe nodesinthehiddenlayerhavebeentakenastheactivationfunctionofalogisticfunction, and in the output layer, the identity function. González-Manteiga et al. Journal of Mathematics in Industry (2018) 8:7 Page 5 of 16 Figure 2 Episode of SO2depicted in one of sampling stations (red) jointly with the prediction provided by the neural network (blue) (Fernández de Castro et al. [6]) The predictor given by the neural network has the following expression: ˆ Xt+6 =o1=L  j=1 ωo 1jfh jθh j+ωh j1Xt–3 +ωh j2Xt with fh j(z)= 1 1+e–z. Theweights {ωh j1,ωh j2,ωo 1j;j=1,...,L}andthetrends{θh j;j=1,...,L}aredeterminedduring the training process, as well as the final Lnumber of hidden layer nodes, that is chosen like the value which neural network provides better results, after having trained networkswithidenticalarchitectureanddifferentvaluesofL.Todesignthetrainingsetofthe neural network it have been considered historical matrices, formerly introduced, suitably adapted. Figure 2shows the forecasts given half an hour before by the neural network with 50 nodes in its hidden layer for an episode depicted in one of the measuring stations. The good behaviour of the forecast (dotted line) can easily be seen. The procedures based on neural networks accurately predict the real one hour mean SO2air quality values (solid line). These models were optimized later with boosting learning techniques (Fernández de Castro and González-Manteiga [4]). 2.1.4 Functional data model TheonehourmeanvaluesofSO2canbetreatedasobservationsofastochasticprocessin continuoustime. Theinterestis,asitwasdiscussedabove,topredictahalf-hourhorizon, so that each of the curves is an interpolated data on half an hour. In this case curves were obtainedbyconsideringsixpentaminutalconsecutiveobservations,withsamplingpoints for each functional data. Therefore, we use random variables with values in Hilbert space H=L2([0,6]) with the form Xt(u)=x(6t+u). González-Manteiga et al. Journal of Mathematics in Industry (2018) 8:7 Page 6 of 16 The following statistical model is considered Xt=ρ(Xt–1)+t,wheretis a Hilbertian strong white noise and ρ:H→His the operator to estimate. For the estimation of ρ, afunctional kernel estimator has been used in the autoregressive Hilbertian of order-one framework.Furthermore,ithasbeenconvenientlyadaptedtheconceptofhistoricalmatrix to the case where the data are curves (Fernández de Castro et al. [5]). 2.1.5 Other approaches designed to predict probabilities Themodelsdescribed,sofar,providepointpredictionsofSO2,butothertechniqueshave also been developed in order to predict probabilities. The aim of these alternative models is to estimate the probability that the series of bihourly SO2measures exceeds a certain level rwith an hour anticipation, namely in our case, we predict P(Zt)=P(Xt+12 >r|Zt) being Zt=(Xt,Xt–Xt–3). To do it additive models with an unknown link function (Roca- Pardiñas et al. [17]) have been used. It has also been considered more complex generalized additive models (GAM) with second-order interaction terms (Roca-Pardiñas et al. [16]). They have shown that the GAMwithinteractionsdetectstheonsetofepisodesearlierthanitdoesGAMonitsown. 2.2 Alternative one-dimensional models: additive models In the statistical literature there is a wide range of one-dimensional models which can be used to predict the levels of SO2. We will focus on the techniques we will use in the next section to construct our multidimensional model: additive models for continuous response. There have been a number of proposals for fitting the additive models. Friedman and Stuetzle [7] introduced a backfitting algorithm and Buja et al. [2] studied its properties. Mammen et al. [12] proposed the so called smooth backfitting by employing projection arguments. Let {(Yt,Zt)}T t=1 be a random sample of a strictly stationary time series, with Ytone-dimensional and Ztq-dimensional following the model: Yt=m(Zt)+t,t∈Z,(1) where {t}is a white noise process and E[t|Zt]=0. Typically,itisassumedthatthefunctionmisadditivewithcomponentfunctionsmj,for j=0,...,q,thus Yt=m0+m1(Z1,t)+···+mq(Zq,t)+t.(2) A generalized kernel nonparametric estimation can be given using smooth backfitting for the functions m1,...,mq(see again the above mentioned papers). Inallthemodelsdescribedaboveitisusuallynecessarytheselectionofaregularization parameter (bandwidth with kernel smoothing, number of neurons in the hidden layer forneuralnetworks,...).Thecalibrationofthisparameterwasdevelopedusingcrossvalidation techniques with the information of the updated Historical Matrix. 3 Methods: multidimensional semiparametric prediction Thenewgoalistoincorporatethepredictionof NOxwithhalfanhourinadvance,aswell as to continuegetting thepredictions of SO2,ashasalreadybeencommented.Theideais González-Manteiga et al. Journal of Mathematics in Industry (2018) 8:7 Page 7 of 16 to generalize the one-dimensional semiparametric approach proposed by García-Jurado et al. [8] taking into account the structure of correlation between the vectorial series that is intended to predict. 3.1 The model Be (Y,Z)=(Yl,Zl), l=0,±1,±2,... avectorialstrictlystationary time series, where Yl is a r-dimensional response series and Zlis a q-dimensional covariables series and, let {(Yt,Zt)}T t=1 be a random sample of (Y,Z). The following model is considered Yt=ϕ(Zt)+Et,(3) where Yt=(Y1,t,...,Yr,t)t,Zt=(Z1,t,...,Zq,t)tand Et=(E1,t,...,Er,t)t. Let us consider two possible structures for the multidimensional residuals series: P1. Each Ek,tis a stationary AR(pk) process of the form Ek,t= pk  i=1 φi kEk,t–i+ξk,tfor all t∈Z,k=1,...,r independentofZt,whereξk,tisawhitenoiseprocesswithvarianceσ2 k,fork=1,...,r. P2. Ethas a VAR(p) structure of the form Et=p  i=1 iEt–i+ξtfor all t∈Z, independent of Zt,wheretheiare fixed (r×r)coefficientsmatricesandξtis a r-dimensional white noise process, i.e. E(ξt)=0,E(ξtξ t)=ξand E(ξtξ s)=0for t=s. Our main objective is to predict Ytusing a sample of size T,κinstants ahead. The prediction of Yt+κis then defined by ˙ Yt+κ=ˆϕκ(Zt)+ ˙ Et+κ,(4) where ˆϕκ(Zt) is a nonparametric estimate of ϕκ(Zt)=E[Yt+κ/Zt]and ˙ Et+κthe prediction given, κinstants ahead, for the residual series constructed as ˆ Et+κ=Yt+κ–ˆϕκ(Zt). 3.2 Estimations We suppose that the model (3) is verified. The first step is to make a nonparametric estimation of ϕindependently foreachof the rcomponents of Yt:ϕ(Zt)=(ϕ1(Zt),...,ϕr(Zt)). Furthermore, we assume that the functions ϕkare additive with component functions ϕj k, for k=1,...,rand j=0,...,q,thus ϕk(Zt)=ϕ0 k+ϕ1 k(Z1,t)+···+ϕq k(Zq,t), k=1,...,r.(5) Therefore, radditive models with qcovariates are estimated using the smooth backfitting technique. We have to take into account that the process Etis not observable since the function ϕis not known. Thus, we have to replace Etby the residuals ˆ Et=Yt–ˆϕ(Zt) González-Manteiga et al. Journal of Mathematics in Industry (2018) 8:7 Page 8 of 16 and use these approximations to Etin the maximum likelihood estimations later defined. To estimate the parametric part of the model, we must consider the two possible error structures proposed above: P1. The parameters φk=(φ1 k,...,φpk k)of the error process {Ek,t}are estimated by standard maximum likelihood methods. In particular, we use a conditional maximum likelihood estimator for every component of the form ˆ φk=argmax φk∈ˆ l(φk), where is a compact parameter space and ˆ lis the conditional log-likelihood given by ˆ lφk,σ2 k=–T 2log(2π)+1 2logσ–2 k–1 2 T  t=pk+1ˆ Ek,t–ˆ Ek,t(φk)/σk2 with ˆ Ek,t(φk)=pk i=1 φi kˆ Ek,t–i. P2. The coefficients matrices (1,...,p)of the r-dimensional error process {Et}are also estimated by generalized maximum likelihood methods (Hamilton [10]). First, we need to establish the following notation: t=[12...p]denote the (r×rp) coefficientsmatrix,letXtbea(rp×1)vectorcontainingplagsofeachoftheelements of Et:Xt t=[Et t–1 Et t–2 ...Et t–p]. Thetheoreticalconditionallog-likelihoodfunctiontobeoptimizedhasthefollow- ing expression: l(,ξ)=–rT 2log(2π)+r 2log–1 ξ–1 2 T  t=1 Et–tXtt–1 ξEt–tXt. Thus the conditional log-likelihood is: ˆ l(ˆ ,ˆ ξ)=–rT 2log(2π)+r 2logˆ –1 ξ–1 2 T  t=1 ˆ Et–ˆ tˆ Xttˆ –1 ξˆ Et–ˆ tˆ Xt. 3.3 Other considerations: the phenomenon of cointegration Sometimesthevectorialprocessescanbecointegrated,soonehastotakeintoaccountthe structure of correlation between the series. The notion of cointegration has been one of the most important concepts in time series since Granger [9] and Engle and Granger [3] that formally developed it. The issue has broad applications in the analysis of economic data as well as several publications in the economic literature. Let Yt=(Y1,t,...,Yr,t)tbe a vector of rtime series integrated of order 1 (I(1)). Ytis said to be cointegrated if a linear combination of them exists that it is stationary (I(0)), i.e., if there exists a vector β=(β1,...,βr)tsuch as βtYt=β1Y1,t+···+βrYr,t∼I(0). The vector βis called the cointegration vector. This vector is not unique since for any scalar cthe linear combination cβtYt=β∗tYt∼I(0). Therefore, normalization is often assumed to identify an unique β. A typical normalization is β=(1,–β2,...,–βr)t. González-Manteiga et al. Journal of Mathematics in Industry (2018) 8:7 Page 9 of 16 Johansen [11] addresses the issue of the cointegration within an error correction model in the framework of vector autoregressive models (VAR). Consider then a general model VAR(p) for the vector of rseries Yt Yt=0Dt+1Yt–1 +···+pYt–p+ξt,t=1,...,T, where Dtcontains deterministic terms (constant, trend, ...). Suppose Ytis I(1) and possibly cointegrated. Then, the VAR representation is not the most suitable representation for analysis because the cointegrating relationships are not explicitly apparent. The cointegrating relationships become apparent if the VAR model is transformed to a vector error correction model of order p(VECM(p)) Yt=0Dt+Yt–1 +1Yt–1 +···+p–1Yt–p+1 +ξt, where =1+···+p–Ir,k=–p j=k+1 j,k=1,...,p–1andYt=Yt–Yt–1.The matrix is called the long-run impact matrix and kare the short-run impact matrices. Moreover,therankofthesingularmatrixprovidesinformationonthenumberofcointegrationrelationsthatexist,i.e.,therankofcointegration.Johansenproposesasequential procedure of likelihood ratio tests to estimate this range. 3.4 Prediction scheme We present now the prediction scheme step by step: 1. Every instant t,ϕκ(Zt)is estimated with the smooth backfitting technique independently for each of rcomponents using the data (Yl,Zl–κ),l=κ+1,...,T. 2. The residuals series ˆ Et+κis computed by ˆ Et+κ=Yt+κ–ˆϕκ(Zt), t=1,...,T–κ. 3. The following step is to make an appropriate adjustment on the model error structure (VECM) and to obtain the prediction κinstants ahead: ˙ ET+κ. 4. The proposed final prediction is given by (4). This scheme is a natural generalization of the one-dimensional prediction models described in Sect. 2.1.1. In the next two sections simulation examples and real data analysis are considered. 4 Results and discussion 4.1 A simulation study To analyze the behavior of the proposed prediction procedure, a simulation study has been performed generating samples from artificial series and making a prediction study to klags using, in all cases, Zt=Yt–1. The following models are considered: Series1. Two independent AR(3) with constant trend: Yt=ϕ+E1,t E2,t, González-Manteiga et al. Journal of Mathematics in Industry (2018) 8:7 Page 16 of 16 12. Mammen E, Linton O, Nielsen J. The existence and asymptotic properties of a backfitting projection algorithm under weak conditions. Ann Stat. 1999;27(5):1443–90. 13. Nadaraya EA. On estimating regression. Theory Probab Appl. 1964;9(1):141–2. 14. Prada-Sánchez J, Febrero-Bande M. Parametric, non-parametric and mixed approaches to prediction of sparsely distributed pollution incidents: a case study. J Chemom. 1997;11(1):13–32. 15. Prada-Sánchez J, Febrero-Bande M, Cotos-Yáñez T, González-Manteiga W, Bermúdez-Cela J, Lucas-Domínguez T. Prediction of SO2 pollution incidents near a power station using partially linear models and an historical matrix of predictor-response vectors. Environmetrics. 2000;11(2):209–25. 16. Roca-Pardiñas J, Cadarso-Suárez C, González-Manteiga W. Testing for interactions in generalized additive models: application to SO2 pollution data. Stat Comput. 2005;15(4):289–99. 17. Roca-Pardiñas J, González-Manteiga W, Febrero-Bande M, Prada-Sánchez J, Cadarso-Suárez C. Predicting binary time series of SO2 using generalized additive models with unknown link function. Environmetrics. 2004;15(7):729–42. 18. Speckman P. Kernel smoothing in partial linear models. J R Stat Soc, Ser B, Stat Methodol. 1988;50:413–36. 19. Watson GS. Smooth regression analysis. Sankhya, Ser A. 1964;26:359–72.