Electricity Spot Prices Forecasting Based on Ensemble Learning

Efficient modeling and forecasting of electricity prices are essential in today’s competitive electricity markets. However, price forecasting is not easy due to the specific features of the electricity price series. This study examines the performance of an ensemble-based technique for forecasting short-term electricity spot prices in the Italian electricity market (IPEX). To this end, the price time series is divided into deterministic and stochastic components. The deterministic component that includes long-term trends, annual and weekly seasonality, and bank holidays, is estimated using semi-parametric techniques. On the other hand, the stochastic component considers the short-term dynamics of the price series and is estimated by time series and various machine learning algorithms. Based on three standard accuracy measures, the results indicate that the ensemble-based model outperforms the others, while the random forest and ARMA are highly competitive.


I. INTRODUCTION
Before the liberalization of the electricity sector, the electrical industry was fully controlled by utility companies, generally state-owned. These utility companies were responsible for the generation, transmission, distribution, and retailing of electricity. In the early 1980s, the global electric industry underwent major changes. The monopolistic electric sector was restructured into a deregulated competitive electricity market. The Public Utility Regulatory Policies Act of 1978 (PURPA) was enacted following the energy crisis of the 1970s to encourage cogeneration and renewable resources and promote competition for electric generation [1]. Soon after, the deregulation process started in Chile in 1982 by introducing the electricity act that dissolved the state-owned monopolistic structure by commercialization and part privatization, followed by large-scale privatization in 1986. In 1990, the British electricity sector started its liberalization process, followed by Norway in 1992. Nowadays, many countries in the world have their own liberalized electricity market.
Liberalization brought important benefits to end consumers such as low prices, more choices, reliable and secure electric supply. However, modeling and forecasting of The associate editor coordinating the review of this manuscript and approving it for publication was Vitor Monteiro . different variables related to electricity market became challenging due to the specific features introduced to this market. As electricity is a commodity that is different from other commodities in terms of its natural and physical qualities [2], electricity price forecasting is very complex due to many characteristics, including non-constant mean, high volatility, spikes or jumps, also known as extreme values [3]. In addition, electricity price prediction is generally required into three time horizons: short, medium, and long-term price prediction [4]. Price prediction for a few hours to a week is generally classified as short-term forecast (STF) which is essential for different market participants. STF helps the producers to plan the electric generation more efficiently. Other market participants require STF to develop superior bidding strategies and maximize profits in day-to-day markets [5]. A medium-term forecast (MTF) refers to a forecast made over a period ranging from a few weeks to several months. The electricity market participants require MTF for a variety of tasks, including maintenance scheduling, generation growth planning, designing investment, fuel contacting, and hedging plans [6]. The long-term price forecast (LTF) generally covers a period ranging from a few months to many years and is commonly used for site planning and selection, investment profitability analysis, and inducing power plant fuel sources [7].
The electricity prices exhibit specific characteristics that include seasonality, volatility, trend, jumps or spikes, etc. [8] and forecasting is more difficult in the presence of these unique characteristics. Due to the highly unpredictable and uncertain nature of the electricity price series, in the past, extensive studies have been made on the problem of electricity prices forecasting using different modeling techniques and procedures. To forecast electricity prices, statistical models such as regression models, time series models, and probabilistic models are commonly used. In the literature, there are a variety of time series models, such as, autoregressive (AR), moving average (MA), autoregressive moving average (ARMA), autoregressive integrated moving average (ARIMA), vector autoregressive (VAR), vector autoregression moving-average (VARMA), seasonal autoregressive integrated moving-average (SARIMA), seasonal autoregressive integrated moving-average with exogenous regressors (SARIMAX), autoregressive conditional heteroskedasticity (ARCH), and generalized autoregressive conditional heteroskedasticity (GARCH) are commonly used for a day-ahead electricity price forecasting [9]- [13]. For next-day electricity prices forecasting, [11] used a nonlinear classical time series model, GARCH, and compared the results to the ARIMA model. They evaluated the proposed methodology on Spanish and California data for hourly electricity prices and observed that the GARCH is significantly more accurate than the ARIMA. One of the most significant steps in power prediction is the shift from a deterministic view to a probabilistic view [14]- [18]. For example, to forecast next-day electricity prices, [17] proposed two modeling strategies, i.e., dynamic regression and transfer function. They used Spanish and Californian market data to evaluate the performance of their proposed methodology. For both markets, the average errors for the weeks analyzed were very low, i.e., approximately 3% and 5%, respectively. Artificial intelligence (AI) models are used to address nonlinear price forecasting issues that are not properly captured by linear time series approaches. The AI or machine learning (ML) approaches have been applied to forecast electricity prices and demand. Artificial neural networks, support vector machine, and long-short term memory (LSTM) models are widely used for electricity price forecasting problems in the literature [19], [20]. For instance, [20] proposed four different deep learning models, i.e., deep neural network (DNN), LSTM, gated recurrent unit (GRU), and convolutional neural network (CNN) for electricity price forecasting. To access the accuracy of these models, they set up a large benchmark study. The results showed that the proposed four models provide the best results with DNN is more accurate among all. In the literature, a combination of different models and algorithms are used to improve forecast accuracy. Multivariate models including functional models are used to forecast oneday-ahead electricity prices [21]- [24]. For instance, [25] used parametric and nonparametric functional models to forecast one-day-ahead electricity prices for five different electricity markets. The results indicated that functional modeling is more efficient in predicting electricity prices than classical univariate and multivariate models. Forecast combination is often described as one of the most effective forecasting strategies. Ensemble modeling is a process in which many individual models are built to predict an outcome, either using a variety of modeling algorithms or different training data set [26]. Nowadays ensemble-based modeling framework is commonly used in the literature to improve forecast accuracy [27]- [32]. For example, [29] proposed an ensemble method which is a combination of empirical mode decomposition(EMD), support vector regression (SVR), and kernel ridge regression (KRR), for forecasting electricity prices in the short-run. Data from the Australian Energy Market Operator (AEMO) is used to access the performance of the EMD-KRR-SVR method. Six benchmark learning algorithms were used to check the accuracy of the proposed method. The results showed that the suggested approach is more efficient and accurate than the benchmark.
Although different models and techniques have been proposed in the literature for electricity price forecasting, none consistently provide accurate forecasting results. In addition, none of these works apply the ensemble approach in the context of component estimation techniques. Furthermore, in general, they do not provide any inferential analysis to test differences in the prediction accuracy of the considered models. Therefore, the main contribution of this research work is the investigation of the ensemble-based approach in the context of component estimation techniques. Furthermore, the forecasting performance is evaluated on a whole year, and the significance analysis of the differences in prediction accuracy is also investigated. In addition, the proposed model can capture the specific features of the electricity price series, leading to higher forecasting accuracy gain.
The rest of the article is organized as follows: Section II describes the methods and models used in this study. The description of the data and an application of the proposed methods and models to the Italian electricity market is given in Section III. The conclusion is given in Section IV.

II. GENERAL MODELING PROCEDURE
The main aim of this research is to forecast one-day-ahead electricity prices by using an ensemble-based modeling structure. As stated, the power price series exhibits a set of different features, and to accurately capture them in a model, the deterministic component filters the price series first, and then the residual component is simulated using a combination of time series and different machine learning algorithms. Once both the components are estimated separately, the final forecasts are obtained by combining the estimates of both deterministic and stochastic components.
Let P d,j denotes the price series for the d th day and j th hour. Then, the price dynamics of P d,j can be described as follows: In other words, the price series P d,j is divided into two major components: ϕ d,j and δ d,j , which represent the VOLUME 9, 2021 deterministic and stochastic components, respectively. The short-run dynamic is accounted by the component δ d,j . On the other hand, the long-run dynamics, different periodicities (yearly and weekly cycles), and calendar influences are all included in the ϕ d,j component. When they are taken into account, ϕ d,j is modeled as follows: where t d,j represents the long-trend, a d,j and w d,j are the annual and weekly cycles, respectively, and c d,j is the calendar effects. Instead of accounting for the daily periodicity in 2, each load period is modeled separately. Finally, the stochastic component δ d,j (residuals) is modeled as which is modeled by different linear and non-linear models.

A. MODELING AND FORECASTING DETERMINISTIC COMPONENT
The estimation and forecasting of the deterministic component ϕ d,j is explained in this section. To capture the different components described in 2, this study considers a semi-parametric approach using the generalized additive modeling technique. That is, penalized smoothing splines are used to model the long-run component (trend) t d,j and the yearly component a d,j , with the series of t d,j and a d,j considered as a functional objects of the time and yearly cycle. Dummy variables are used to describe weekly cycle w d,j and bank holidays c d,j . In the case of weekly periodicity, seven dummy variables are required, i.e., w d,j = 7 i=1 γ i I d,j where I d,j = 1 if d is the i th day of the week, and 0 otherwise. In the case of bank holidays, only two dummy variables are used, i.e., c d,j = 2 k=1 α k I d,j , where I d,j = 1 if d indicates a bank holiday and 0 otherwise. The parameters γ i and α k are estimated using the ordinary least squares (OLS) methods. Once these components are estimated, their one-day-ahead forecasts are straightforward as all these components are the deterministic functions of time or calendar conditions.

B. MODELING AND FORECASTING STOCHASTIC COMPONENT
Following the estimation of the deterministic components using a semi-parametric approach, the stochastic (residual) component δ d,j is modeled using a combination of autoregressive moving average (ARMA) and different machine learning models, such as neural network auto-regressive (NAR), random forest (RF), support vector regression (SVR), gradient boosting machine (GBM), and an ensemble model based on these methods. Once these models are estimated and a oneday-ahead forecast is obtained for each model, the final oneday-ahead out-of-sample forecast is obtained by combining both the deterministic and stochastic components as follows.
The details about each model used for the estimation of the stochastic component is given in the following.

1) AUTOREGRESSIVE MOVING AVERAGE (ARMA) MODELS
A powerful tool for modeling univariate time series is an autoregressive moving average (ARMA) model. In general, an ARMA model is a combination of autoregressive (AR) and moving average (MA) models. The mathematical form for an ARMA of order (s, r), abbreviated as ARMA(s, r), is as follows.
where β represents the constant term (intercept), ψ j (j = 1, 2, · · · , s) and ϑ i (i = 0, 1, 2, · · · , r) are the AR and MA parameters, respectively, and ζ d−r,j is a white noise series with zero mean and variance σ 2 ζ . In general, different information criteria or auto-correlation function (ACF) and partial auto-correlation function (PACF) are used to specify the order of s and r.

2) ARTIFICIAL NEURAL NETWORK AUTOREGRESSIVE (NNAR)
Generally, a neural network is a network or circuit of neurons. An artificial neural network is made up of artificial neurons or nodes. Modeling complex nonlinear relationships between the response variable and its predictors are possible using neural network models. A feedback neural network is assembled with delayed time series values as input and a hidden layer with dimension nodes. There are at least three layers of nodes in an autoregressive neural network (NNAR): an input layer, a hidden layer, and an output layer. The outputs of one layer are used as inputs to the next. The nonlinear autoregressive neural network can be trained to forecast a time series δ d,j given its previous values δ d−1,j , δ d−2,j , · · · , δ d−p,j , which are referred as feedback delays, with d is the time delay parameter. The term NNAR(p,k) implies that the hidden layer contains p delayed inputs and k nodes. An NNAR(p,0) is the same as an ARMA(p,0), but it does not include the limitations of the parameters that ensure stationarity. A nonlinear autoregressive neural network used for time series prediction can be described as follows.
Since the γ (.) function for the next day is unknown, the neural network training aims to approximate it by optimizing the network weights and neural bias. Consequently, the following equation precisely specifies an NNAR model.
where q is the number of entries, p is the number of hidden layers with activation function θ, and φ k i is the parameter corresponding to the weight of the connection between the input unit k and the hidden unit i, α i is the weight of the connection between the hidden unit i and the output unit, and φ i and α are the constants that correspond to the hidden unit i and the output unit, respectively. In this work, NNAR(2,3) is used which implies that the hidden layer contains 2 delayed inputs and 3 nodes.

3) RANDOM FOREST (RF)
Random decision forests, often known as random forests, are ensemble learning methods for classification and regression that work by training a large number of decision trees. Breiman et al. [33] presented the procedures of random forest (RF), classification and regression tree (CART). In 1996, Breiman proposed another significant RF technique called Bagging [34]. The RF is a well-known machine learning algorithm that falls within the supervised learning category. The RF algorithm is applied in a wide range of areas, including finance, stock trading, health care, and e-commerce. It produces a forest out of a collection of decision trees that are often trained using the bagging method. The RF algorithm determines the outcome based on the decision trees' predictions. It forecasts by averaging the output of various trees. The precision of the result improves as the number of trees grows. The higher the number of trees in the forest, the more accurate it is and the problem of over-fitting is avoided. The random forest training algorithm uses the common technique of bootstrap aggregation, or bagging, to train tree learners. Bagging repeatedly (M times) takes a random sample with replacement of the training set and fits trees to these samples given a training set Y = y 1 , . . . , y n with responses X = x 1 , . . . , x n . For m = 1, . . . , M , the algorithm works as follows.
1. Sample m training data points from X and Y with replacement and denote them Y m and X m .
2.Train a classification or regression tree f m on Y m and X m . After training, summing the predictions from all the various regression trees on y can be used to make predictions for unseen samples y .f By reducing the variance of the model without increasing the bias, this bootstrap approach improves the performance of the model. In the case of classification trees, the majority vote is used.

4) GRADIENT BOOST MACHINES (GBM)
The gradient boost machine (GBM) is a type of ensemble learning that uses a sequential learning process to create a good classification or regression model. Initially, the data are fitted to a regression tree, and predictions and the initial residual are produced based on this information. A new model is fitted to the previous residual, followed by a new forecast, which includes the initial forecast and finally a new residual. This procedure is continued iteratively until a convergence condition is satisfied. At each iteration, a new model is fitted to the data to compensate for the deficiencies of the prior model. VOLUME 9, 2021 The GBM algorithm works best when the contribution of the additional decision tree is minimized at each step of the iteration by using shrinkage parameter β, also known as the learning rate. In the GBM, the idea behind the shrinkage procedure is that more small passes provide greater precision than fewer large passes. The learning parameter β has a range of values between 0 and 1; the smaller the value, the more accurate the model. However, because the value of β is inversely proportional to the number of iterations, choosing a stronger shrinkage (smaller β) implies a higher number of iterations to achieve convergence. The following steps provide a simplified illustration of this algorithm: 1. Sets the hyper parameters of GBM: depth of the decision trees τ , the number of iterations p, the shrinkage parameter β, and the fraction of the subsample, n.
2. Set ζ 0 = x and f = 0 for the residual. The mean value of x has also been proposed as an initial estimate of f .
3. For p = 1,2,. . . ,P, follow: a. Select a subsample {xi, yi} N * at random from the entire training dataset, where N * is the number of data points corresponding to the fraction n.
b. Fit a decision tree f p of depth τ to the residual ζ p−1 value using {xi, yi} N * .
c. Add the decision tree to the model to update f .
d. Update the residual There are four hyperparameters in the GBM: (1) decision tree depth (τ = 1), which also affects the model's maximum interaction order, (2) number of iterations (p), which is same as the number of decision trees, (3) learning rate (β = 0.001), which is normally a small positive value between 0 and 1, (4) fraction of the subsample (n), which is the percentage of data used in each iteration step.

5) SUPPORT VECTOR REGRESSION (SVR)
Support vector machine can also be applied as a regression approach while retaining all of the algorithm's key properties (maximal margin) [35]. With a few minor exceptions, the support vector regression (SVR) uses the same classification concepts as the SVM. The SVR involves finding support vectors adjacent to a hyperplane that optimizes the gap between two-point classes based on the difference between the target and a threshold. The SVR uses kernel functions to calculate the similarity between two observations to handle non-linear situations. The capability of the SVR to capture predictor non-linearity and then use it to improve forecasting situations is one of the key advantages of using it. The SVR helps to determine how much error is acceptable in the model. Support vector regression equation with kernel function can be written as The kernel function, K (y i , y), represents the inner product, while α is accommodated within the kernel function.
where y i − y k 2 is the Euclidean distance between the two feature vectors squared.

6) ENSEMBLE MODEL
Ensemble methods are strategies that combine multiple models to achieve better results. In most cases, ensemble approaches provide more accurate results than a single model. In many machine learning competitions, the winning solutions used ensemble approaches. In this work, five models are used to make an ensemble namely, autoregressive moving average (ARMA), neural network auto-regressive (NAR), random forest (RF), support vector regression (SVR), and gradient boosting machine (GBM). To make an ensemble model weighting scheme is applied. A weighted average ensemble is a method for allowing different models to contribute to a forecast in proportion to their level of confidence or estimated performance. The contribution of each member to the final forecast is weighted by the model's performance in a weighted ensemble, which is an extension of a model averaging ensemble. The model weights are small positive values and the sum of all weights equals one, indicating the percentage of trust or expected performance from each model.

7) NAÏVE MODEL
This study used the benchmark model proposed by [36]. Although it is a naïve model, the forecasting accuracy of this approach is much better when compared with other naïve models. This method forecasts a day based on the information contained in all the preceding days. For example, if we are interested in predicting a Saturday, we select the day before Saturday, i.e., Friday. We then compare this Friday with all available Fridays in the training set and calculate their errors. Based on the error, we decide which Friday is more similar to the current Friday. Once a similar Friday is identified, we select the next day to the selected Friday, i.e., Saturday, and use it as a forecast for the Saturday we are interested in. We do the same procedure for all days of the week and for the whole year.

III. OUT-OF-SAMPLE PRICE FORECASTING
The electricity price data from the Italian electricity market (IPEX) is used in this work to evaluate the forecasting performance of different models. The IPEX is an important electricity market that is fully liberalized since 2007. For this study, the data includes 52584 observations over 2191 days, ranges from January 1st, 2014 to December 31st, 2019. The data collection contains 24 observations for each day corresponds to a load period and the prices are given in Euro/MWh. The descriptive statistics for the data under study are given in Table 1. The table provides the minimum, the maximum, and the mean prices as well as the standard deviation of the weekdays' prices. The table shows that the minimum price reaches 0.10 for Wednesday and Sunday while the maximum price is 324.20 for Tuesday. The mean prices for the weekdays range from 48.70 to 59.34. The last column provides the standard deviation values for each day of the week, while the last row describes the descriptive statistics for the whole electricity price time series. For modeling and forecasting, the data set is divided into two groups: For model identification and estimation, we used data from January 1st, 2014 to December 31st, 2018 (43824 hourly observations/ 1826 days), and January 1st, 2019 to December 31st, 2019 (8760 hourly observations/ 365 days) for analyzing the models' one-day-ahead out-of-sample forecasting accuracy using the expending window technique. The electricity price time series is plotted in Figure 2 where the vertical red dotted line divided the estimation and out-of-sample forecasting periods. The main reason for using a whole year as an outof-sample forecasting period is to assess the performance of each model in every type of situation during a year. Since the prices vary throughout the year, with some periods producing extreme prices, evaluating models' performance on a complete year is a more realistic approach. The deterministic component is first estimated and forecasted to obtain the stochastic (residual) part. For the stochastic component, all the models are used to estimate and forecast the stochastic component. This study considers three standard accuracy measures, namely mean absolute error (MAE), root mean squared error (RMSE), and Pearson correlation (R), to assess the model's forecasting performance. Mathematically, they are defined as The observed and forecasted prices for the dth day (d = 1, 2, . . . , 365) and the j th load period (j = 1, 2, . . . , 24) are P d,j andP d,j , respectively. the MAE of the ensemble model. Moreover, the forecasting accuracy of the naïve model is much lower than all other models used in the study. This shows that the current modeling framework is efficient in predicting one-day-ahead electricity prices.
Several methods are available in the literature to compare the performance of different forecasting models [37]. In this work, the Diebold and Mariano (DM) [38] test is used to determine the significance of the differences between the results reported in Table 2, and the results are listed in Table 3. The null hypothesis of the DM test corresponds to equal forecast accuracy, while the alternative hypothesis states that the model in the row is more accurate than in the column. The DM results listed in Table 3 indicate that the proposed ensemble model is statistically significant than all other models. In addition, the proposed modeling structure is also statistically significant than the naïve model.
Finally, the obtained results are compared with the results provided in the literature and found to be highly comparative. It is worth mentioning that such comparisons are very difficult as different authors use different indicators, forecasting periods, and forecasting horizons. For example, the work of [39] reported an MAE of 8.58, whereas we obtained a value of 5.05 with the ensemble model. Ref. [40] reported an RMSE of 11.58 obtained with an ARX-EGARCH model, whereas we obtained 6.66 with the ensemble model. Ref. [41] obtained an RMSE of 16.72 and 15.79 using ARMA and GARCH models, respectively, which is considerably higher than our value of 6.66.

IV. CONCLUSION
The forecasting problem of electricity prices in the recently liberalized market is analyzed in depth in this work. In competitive electricity markets, efficient modeling and forecasting of electricity prices are essential as forecasts are necessary for risk management, trading, and future planning because prices and demand are determined one day before physical delivery. However, price forecasting is not easy as electricity prices can have complex characteristics, such as high volatility, seasonality, calendar impact, non-stationarity, nonlinearity, mean reversion, etc. Thus, this research work aims to propose an ensemble-based technique to improve prediction accuracy by combining time series and different machine learning algorithms. To this end, the deterministic component filters the price series first, and then the residual component is modeled using an ensemble-based method to accurately capture the specific properties of electricity prices in the model. The deterministic component that includes long-term trends, annual and weekly seasonality, and bank holidays, is estimated using a semi-parametric approach. The stochastic component, on the other hand, considers the short-term dynamics of the price series and is modeled by the autoregressive moving average (ARMA), and various machine learning algorithms, including neural network autoregressive (NNAR), random forest (RF), support vector regression (SVR), gradient boosting machine (GBM). After modeling the deterministic and stochastic components individually, the final prediction is obtained by combining the estimates of both components. The proposed modeling approach is presented together with an application to the Italian electricity market (IPEX). To evaluate the performance of our models, data from the period January 1st, 2014 to December 31st, 2019 is used, whereas from January 1st, 2014 to December 31st, 2018 is used for model estimation while the entire 2019 year is for the one-day-ahead out-of-sample forecast. The mean absolute error (MAE), root mean square error (RMSE), and Pearson's correlation coefficient (r) are used to assess the effectiveness of the ensemble-based methodology.
The results suggest that the proposed ensemble method is efficient in predicting electricity prices. The ensemble model produced the MAE value of 5.05, whereas the ARMA and RF produced 5.18 and 5.23, respectively. The results of the NNAR with MAE 5.29, GBM with MAE 5.30, and SVR with MAE 5.39 are also good, but not satisfactory as compared to the MAE of the ensemble model. Moreover, the forecasting accuracy of the nave model with the MAE of 6.66 is much lower than all other models used in the study. The DM test validated the significance of the ensemble model compared to the rest. Lastly, as the ensemble-based model outperforms the others, the RF and ARMA produce better results than other models used in this study. In conclusion, the proposed modeling framework is statistically significant in forecasting electricity prices. As the current study considers only the IPEX, in the future, it can be extended to other markets to assess the performance of the proposed models.
NADEELA BIBI is currently pursuing the M.Phil. degree in statistics with the Department of Statistics, Quaid-i-Azam University (QAU), Islamabad, Pakistan. Her research interests include time series analysis, forecasting electricity market variables, and applied statistics.
ISMAIL SHAH received the master's degree from Lund University, Sweden, and the Ph.D. degree from the University of Padua, Italy. He is currently working as an Assistant Professor with the Department of Statistics, Quaid-i-Azam University, Islamabad, Pakistan. His research interests include functional data analysis, time series analysis, regression analysis, energy economics, and applied and industrial statistics. He is also working as an Editor of the Journal of Quantitative Methods.
ABDELAZIZ ALSUBIE received the Ph.D. degree in statistics from North Dakota State University, in 2020, under the supervision of Prof. Rhonda Magel. He is currently working as an Assistant Professor of statistics with the Department of Basic Sciences, Saudi Electronic University. His research interests include nonparametric, inference under order restriction, and simulation studies.
SAJID ALI received the Ph.D. degree in statistics from Bocconi University, Milan, Italy. He is currently an Assistant Professor with the Department of Statistics, Quaid-i-Azam University (QAU), Islamabad, Pakistan. His research interests include time series analysis, Bayesian inference, construction of new flexible probability distributions, and process monitoring.
SHOWKAT AHMAD LONE received the M.Sc. degree in statistics from the University of Kashmir, in 2013, and the Ph.D. degree in statistics from Aligarh Muslim University, India, in 2017. For more than seven years, he has been teaching courses in statistics, mathematics, management science, economics, and other quantitative models. He is currently an Assistant Professor of mathematics and statistics with Saudi Electronic University (SEU), Saudi Arabia. He has authored textbook chapters in reliability engineering and statistics. His research interests include the area of reliability engineering, computer science, and related fields. He has published research articles in reputed international journals in the fields of mathematical and statistical sciences. He is a referee for several international journals in the field of pure and applied statistics.