Forecasting the Power Generation Mix in Italy Based on Grey Markov Models

1. Introduction

Grey Systems Theory (GST) was initiated by Julong Deng in 1982 [1] and further developed in [2]. It represents an approach aimed at forecasting systems characterized by small sample sizes and poor information. This methodology excels in managing systems where only fragments of information are known by identifying, unearthing, and refining the valuable data present. Grey systems theory has proven to be an effective strategy for handling uncertainty in information. In applied sciences, it’s common to encounter systems that are not fully understood and are supported by minimal data sets and vague information, highlighting the need to apply this theory.

Since its inception, GST has witnessed consistent expansion in terms of cross-disciplinary applications and in the diversity of its methodologies. In particular, the proposal of novel approaches and techniques for addressing data that are unknown, incomplete, or inadequate is at the core of recent developments.GST has been employed across a broad spectrum of fields, for example, in predictive scientometrics [3], energy and environmental emissions [4,5], materials science [6], engineering disciplines [7], the development of electric vehicles [8], environmental studies [9], management of water leaks [10], within the fields of economics and social sciences [11], and additional examples are illustrated in [12,13].

Generally, the Grey model $G(M,N)$ comprises M-order partial differential equations and N variables. The principle of grey modelling is founded on the assumption that a process’s underlying mechanisms can be abstracted by a specific, standard partial differential equation, which mirrors the method’s degree of complexity. The estimation of model parameters is conducted through the analysis of accident data. This modelling framework is distinguished by its ability to use just four pieces of data to project future events [14]. In the literature, two single-variable grey models are predominantly utilized: the first-order $G(1,1)$ and the second-order $G(2,1)$.

In practice, observed data are frequently subject to various random influences, leading to data series that exhibit unpredictability and volatility. The $GM(1,1)$ forecasting model is an effective tool for data sequences with minimal data, but it struggles with sequences that exhibit significant fluctuations [15]. To overcome this issue and to increase forecast precision, the Grey Markov Chain Model $MCGM(1,1)$ was created to correct the residual errors produced by $GM(1,1)$. This model merges the first-order single-variable $GM(1,1)$ grey forecasting model with a Markov-chain ($MC$) model that adjusts for residuals [16,17]. This integration of models has been widely applied to datasets that are small in size but characterized by randomness and fluctuation [18,19,20]. The core concept behind the $MCGM(1,1)$ model involves identifying discrete states for the residual errors from $GM(1,1)$ and computing Markov probability transition matrices to grasp the statistical properties of these random residuals. Subsequently, adjustments are made based on these Markov matrices to refine the forecasting accuracy of the original $GM(1,1)$ model.

In this article, we have chosen to retrace the steps of the model described above for forecasting the energy mix of primary energy consumption for Italy up to 2032.

The importance of this topic lies in the strategic planning and policy-making to meet environmental goals and energy security. Forecasting the energy mix enables the assessment of future energy scenarios, including the adoption of renewable energy sources, and helps in evaluating their impact on carbon emissions reduction and energy security. A recent application study to the Australian framework is discussed in [21].

For Italy, as part of the European Union (EU), forecasting the energy mix is crucial for aligning with the EU’s ambitious targets for reducing greenhouse gas emissions and increasing the share of renewable energy in total energy consumption by 2030 and beyond. Understanding the future energy mix aids in planning for the necessary infrastructure developments, such as grid upgrades and energy storage solutions, to accommodate an increased share of variable renewable energy sources like wind and solar power, see e.g. [22].

Moreover, forecasting the energy mix is essential for identifying potential challenges and opportunities within the energy transition. For instance, it can highlight the need for Italy to diversify its energy sources to reduce dependence on imported fossil fuels, thus enhancing energy security. It also allows for the identification of economic opportunities related to the energy transition, including job creation in the renewable energy sector and potential for technological innovation.

Furthermore, an accurate forecast of the energy mix is vital for ensuring the reliability and stability of the energy supply. As Italy moves towards a more sustainable energy system, understanding the balance between supply and demand, the integration of renewable energy and the role of energy efficiency becomes increasingly important. This insight is critical for maintaining a stable energy supply while meeting environmental targets.

In conclusion, forecasting the energy mix is of paramount importance for Italy as it navigates its energy transition. It supports strategic planning, policy development, and infrastructure investment decisions aimed at achieving a sustainable, secure, and competitive energy system.

Given the absence of prior applications of this model for this specific purpose in the literature, we have endeavoured to outline its implementation steps towards the applied results.

This study provides an insightful look into Italy’s future energy composition, using a forecasting model that integrates the Grey system theory with Markov chains to address the uncertainties in energy prediction. The obtained results reveal a significant trend towards the adoption of renewable energy sources while indicating that fossil fuels will still play a dominant role in Italy’s energy mix by 2030. This dual dependence underscores the complexity of transitioning to a sustainable energy framework and highlights the critical need for Italy to enhance its efforts in renewable energy development. The results serve as a call to action for policymakers, suggesting that while strides are being made towards a greener energy future, more aggressive policies and investments are required to reduce fossil fuel dependency significantly. This study not only sheds light on the potential energy landscape of Italy in the coming years but also emphasizes the importance of innovative forecasting techniques in planning and policy-making for energy transition.

The paper is organised as follows. In Section 2, we present the model. Specifically, in the SubSection 2.1, we illustrate the steps of the Grey Model GM(1,1), and in the SubSection 2.2, we explain how to apply the Markov chain for error adjustment. Subsequently, in Section 3, we introduce our case study focused on forecasting the consumption of primary energy from each energy source for Italy until 2032 and in Section 4 the conclusions of our work are given.

2. The Model

The GM(1,1) model employs a first-order differential equation to describe the forecasting target, and it stands as a primary tool in grey system forecasting. Additionally, the enhanced MCGM(1,1) approach, merging the principles of the GM(1,1) model with those of the Markov chain, delivers improved statistical assessments for datasets that demonstrate significant fluctuations.

2.1. GM(1,1) Model

The GM(1,1) model is widely recognized for its ability to generate reliable short-term forecasts using a limited dataset of non-negative values. In this section, we will discuss the steps required to implement this model.

Let $X^{\left(0\right)}$ be the vector representing the original data sequence containing n observations:

$${\mathit{X}}^{\left(0\right)}={(x_1^{\left(0\right)},x_2^{\left(0\right)},...,x_n^{\left(0\right)})}^T.$$

(1)

Firstly, we perform a one-time accumulated generating operation in the following way

$$x_k^{\left(1\right)}=\sum _{i=1}^k{x}_i^{\left(0\right)},k=1,2,\dots ,n,$$

(2)

in order to create the ${\mathit{X}}^{\left(1\right)}$ array

$${\mathit{X}}^{\left(1\right)}={(x_1^{\left(1\right)},x_2^{\left(1\right)},...,x_k^{\left(1\right)},...,x_n^{\left(1\right)})}^T.$$

(3)

It should be observed that the first element of vector in Eq. (1) corresponds to the first element of vector in Eq.(2), symbolically represented as $x_1^{\left(0\right)}=x_1^{\left(1\right)}$.

Then, set $Z^{\left(1\right)}=(z_2^{\left(1\right)},z_3^{\left(1\right)},...,z_n^{\left(1\right)})$ be the vector of background values computed as follows:

$$z_k^{\left(1\right)}=\frac{x_k^{\left(1\right)}+x_{k-1}^{\left(1\right)}}{2},k=2,\dots ,n.$$

(4)

This is equivalent to stating that the general element $z_k^{\left(1\right)}$ is the arithmetic mean of the two neighbouring data points.

The $GM(1,1)$ equation is defined as follows

$$x_k^{\left(0\right)}+a{z}_k^{\left(1\right)}=b$$

(5)

where a is the developing coefficient and b is the Grey effect.

Next, the cumulative values $x_1^{\left(1\right)},x_2^{\left(1\right)},\dots ,x_n^{\left(1\right)}$ can be approximated by the first-order Grey differential equation

$$\frac{d{X}^{\left(1\right)}}{dt}+a{X}^{\left(1\right)}=b,$$

(6)

Based on the principles of Least Square Estimation (LSE), the estimation of parameters $\widehat{a}$ and $\widehat{b}$ can proceed in the following manner:

$${[\widehat{a},\widehat{b}]}^T={\left({\mathit{B}}^{\mathit{T}}\mathit{B}\right)}^{-1}{\mathit{B}}^{\mathit{T}}\mathit{y},$$

(7)

where

$$\mathit{B}=\left[\begin{array}{cc}-z_2^{\left(1\right)}& 1\\ -z_3^{\left(1\right)}& 1\\ ⋮& ⋮\\ -z_n^{\left(1\right)}& 1\end{array}\right],\mathit{y}=\left[\begin{array}{c}{x}_2^{\left(0\right)}\\ x_3^{\left(0\right)}\\ ⋮\\ x_n^{\left(0\right)}\end{array}\right].$$

(8)

To calculate an estimate of the Standard Error $\left(SE\right)$ associated with the parameters $\widehat{a}$ and $\widehat{b}$ we compute

$$S{E}^2=\widehat{Var}\left[\begin{array}{c}\widehat{a}\\ \widehat{b}\end{array}\right]=\left(\frac{{\sum }_{k=1}^n{\left(x_k^{\left(0\right)}-{\widehat{x}}_k^{\left(0\right)}\right)}^2}{n-2}\right){\left({\mathit{B}}^T\mathit{B}\right)}^{-1},$$

(9)

Finally, the solution of [6] led us to the time response function

$${\widehat{x}}_k^{\left(1\right)}=\left(x_1^{\left(0\right)}-\frac{\widehat{b}}{\widehat{a}}\right)e^{-\widehat{a}(k-1)}+\frac{\widehat{b}}{\widehat{a}},k=1,2,\dots$$

(10)

Therefore, the predicted values of $X^{\left(0\right)}$ are:

$${\widehat{x}}_k^{\left(0\right)}={\widehat{x}}_k^{\left(1\right)}-{\widehat{x}}_{k-1}^{\left(1\right)}=\left(1-e^{\widehat{a}}\right)\left(x_1^{\left(0\right)}-\frac{\widehat{b}}{\widehat{a}}\right)e^{-\widehat{a}(k-1)},k=2,\dots ,$$

(11)

given the assumption that the first value is known, that is ${\widehat{x}}_1^{\left(0\right)}=x_1^{\left(0\right)}$.

2.2. Markov Chain Residual Modification

Observed data often fluctuate due to various random influences. Consequently, accurate forecasts with the GM(1,1) model are difficult to obtain. To overcome this problem, we integrate a Markov-chain approach to refine the residual errors. In particular, the MCGM(1,1) model merges the grey forecasting approach of GM(1,1) with a Markov-chain mechanism to adjust for residual errors. The core concept behind the MCGM(1,1) involves first identifying discrete states for the residual errors from GM(1,1), then calculating Markov probability transition matrices to understand the statistical behaviour of these residuals, and finally applying corrections based on these matrices. This process enhances the original forecasting accuracy of the GM(1,1) model. The step-by-step residual modification is described in the current section.

First, we define which residuals to be used. There are examples in the literature that show the use of absolute errors [23] and others that show the use of relative errors [20,24]. In our investigation, we focus on relative errors represented in the form

$$r{\epsilon }_k=\frac{{\widehat{x}}_k^{\left(0\right)}-x_k^{\left(0\right)}}{x_k^{\left(0\right)}},k=1,2,\dots ,n.$$

(12)

Let ${\epsilon }_k=({\epsilon }_1,{\epsilon }_2,...,{\epsilon }_n)$ be the chosen series of residual values, ${\epsilon }_{min}$ be the minimum and ${\epsilon }_{max}$ be the maximum. We split the range of residuals $[{\epsilon }_{min},{\epsilon }_{max}]$ into q intervals of the same length. The generic interval is denoted by

$$r_w=\left[{\epsilon }_{min}+(w-1)\frac{{\epsilon }_{max}-{\epsilon }_{min}}{q};{\epsilon }_{min}+w\frac{{\epsilon }_{max}-{\epsilon }_{min}}{q}\right],w=1,\dots ,q.$$

(13)

Each of them is assumed to be a state. In particular, state 1 is the interval with the lower bound equal to ${\epsilon }_{min}$, and state q is the one with ${\epsilon }_{max}$ as the upper bound.

This means that the sequence of residual values is converted into a sequence of states assuming values in the set $E=\{r_1,r_2,\dots ,r_q\}$. The state of residual ${\epsilon }_k$ depends on which interval belongs to. The specific number of state divisions is not rigidly defined; rather, it is determined based on comprehensive factors such as the sample size and the error range of the fit [25]. Classifying them into 3 to 5 states is typically deemed suitable [26].

Now, take $c_w$ as the notation for the representative value of each state $r_w$, $w=1,2,\dots ,q$ whose lower bound is $l_w$ and upper bound is $u_w$. In some circumstances, simple choices are considered. For example, the formula proposed in [23] is a weighted average of the lower and upper bounds of the intervals:

$$c_w={\alpha }_w{l}_w+(1-{\alpha }_w)u_w,w=1,2,\dots ,q,$$

(14)

where ${\alpha }_w$ are the weights that can be assigned to the extremes of each interval. In general, various choices are possible by producing statistical estimations of the representative values of the states of the Markov process. This is commonly done in diverse application settings, see for example [27,28]. Accordingly, we consider the following estimation formula:

$$c_w=\frac{{\sum }_{k=1}^n{\mathbb{1}}_{\{{\epsilon }_k\in r_w\}}·\left({\epsilon }_k\right)}{{\sum }_{k=1}^n{\mathbb{1}}_{\{{\epsilon }_k\in r_w\}}}.$$

(15)

The next assumption is to consider this series of states as generated by a Markov chain ${\left\{R_k\right\}}_{k\in \mathbb{N}}$. Hence, for every time $k\in \mathbb{N}$, it results that

$$\mathbb{P}[R_{k+1}=c_j|R_k=c_i,R_{k-1}=c_h,...]=\mathbb{P}[R_{k+1}=c_j|R_k=c_i]=:p_{ij}\equiv p_{ij}^{\left(1\right)}.$$

(16)

The above probabilities define the so-called one-step transition probability matrix, $\mathbf{P}={\left(p_{ij}\right)}_{i,j\in E}$.

Now, we introduce the transition probability matrix ${\mathit{P}}^{\left(m\right)}$ as an $q\times q$ matrix with elements $p_{ij}^{\left(m\right)}$, where i and j denote the rows and columns respectively, and m counts the number of the steps. Indeed, ${\mathit{P}}^{\left(m\right)}$ is called the m-steps transition probability matrix and represents the m-th power of the matrix $\mathit{P}$.

$${\mathit{P}}^{\left(m\right)}=\left[\begin{array}{cccc}{p}_{11}^{\left(m\right)}& p_{12}^{\left(m\right)}& \cdots & p_{1q}^{\left(m\right)}\\ p_{21}^{\left(m\right)}& p_{22}^{\left(m\right)}& \cdots & p_{2q}^{\left(m\right)}\\ ⋮& ⋮& \ddots & ⋮\\ p_{q1}^{\left(m\right)}& p_{q2}^{\left(m\right)}& \cdots & p_{qq}^{\left(m\right)}\end{array}\right]$$

(17)

According to the Markov chain assumption, it results that

$$p_{ij}^{\left(m\right)}=\mathbb{P}[R_{k+m}=c_j|R_k=c_i].$$

(18)

Each element $p_{ij}^{\left(m\right)}$ represents the probability of $R_{k+m}=c_j$ given that $R_k=c_i$, i.e. the residual process transits from state $c_i$ to state $c_j$ in m steps. The one-step transition probability matrix is estimated as follows:

$$p_{ij}^{\left(1\right)}=\frac{t_{ij}^{\left(1\right)}}{t_i},i,j=1,2,\dots ,q.$$

(19)

The term $t_{ij}^{\left(1\right)}$ represents the total number of transitions from the i-th to the j-th state after $m=1$ step, and $t_i$ is the count of occurrences of state i within the residual sequence. To obtain the m-step transition matrix, the one-step matrix is raised to the m-th power.

Each row in ${\mathit{P}}^{\left(m\right)}$ has elements that sum to one. If the sum of a row’s elements for state i is nil because we never observed a residual belonging to the i-th interval $r_i$, one choice could be to directly assign $p_{ii}^{\left(1\right)}=1$. This choice implies that state $r_i$ functions as an absorbing state.

Due to limited data, another problem may arise in the application of this model. It may be possible that a state of the process is only visited once and at the last observation of the time series of the residuals. In this eventuality, we are unable to estimate the transition probabilities over the row of the transition matrix which corresponds to the considered state. Here, the choice to characterize that state as an absorbing one poses a major risk when forecasting over a time span exceeding the observation time. In such a scenario, rendering state $r_i$ absorbing exclusively dictates the predictions, thereby introducing distortions. This is because it is unrealistic to assume that the process remains in the residual state $r_i$ until the end of the predictive interval.

Therefore, our recommendation is to assign to each element over the row an equal probability. This equates to the random initialization of errors within the system, subsequently allowing them to develop following the observed frequencies via the Markovian framework.

For example, this is equivalent to saying that for a number of states equal to $r=4$, each row element would be equal to $0.25$. In this way, when the chain reaches this particular state, it will not be absorbed by the starting state and it will be equally likely to reach in the next visit any other state.

The Markov chain method is employed to adjust the initial values obtained from the GM(1,1) model predictions revising it by adding the effect of a predicted residual. We can now describe this procedure, as follows.

Let us assume that at the current time k the residual process is in the state i, that is ${\epsilon }_k\in r_i=[l_i,u_i]$. The GM(1,1) model produces forecasting ${\widehat{x}}_{k+m}^{\left(0\right)}$ of the variable $x^{\left(0\right)}$ for future times $k+m,m=0,1,2\dots$. These forecasts are adjusted according to the formula

$${\widehat{\widehat{x}}}_{k+m}^{\left(0\right)}={\widehat{x}}_{k+m}^{\left(0\right)}(1-{\widehat{\epsilon }}_k^{\left(m\right)}),k=1,2,\cdots ,n.$$

(20)

The prediction ${\widehat{\epsilon }}_k^{\left(m\right)}$ can be computed conditionally on the state of the residual process $R_k$. Thus, if the residual at time k is in state $r_i$, the prediction is given by:

$${\widehat{\epsilon }}_k^{\left(m\right)}=p_{i1}^{\left(m\right)}{c}_1+p_{i2}^{\left(m\right)}{c}_2+\cdots +p_{ir}^{\left(m\right)}{c}_r,k=1,2,\cdots ,n.$$

(21)

The prediction formula (20) expresses a m-steps forecasting.

Thus, for $m=0$, we obtain a contemporary adjustment according to the state $r_i$ given by the relation

$${\widehat{\widehat{x}}}_k^{\left(0\right)}={\widehat{x}}_k^{\left(0\right)}(1-c_i),k=1,2,\cdots ,n.$$

(22)

which is obtained setting $m=0$ into (20) and (21). Hence, the Markov chain corrects the initial prediction according to the representative value of the residual state at time k.

For $m=1$ and for the residual at time k being at state $r_i$, we obtain

$${\widehat{\widehat{x}}}_{k+1}^{\left(0\right)}={\widehat{x}}_{k+1}^{\left(0\right)}(1-{\widehat{\epsilon }}_k^{\left(1\right)})={\widehat{x}}_{k+1}^{\left(0\right)}(1-\sum _{j=1}^q{p}_{ij}{c}_j),k=1,2,\cdots ,n.$$

(23)

In this case, the GM(1,1) forecast ${\widehat{x}}_{k+1}^{\left(0\right)}$ is corrected according to the factor $1-{\sum }_{j=1}^q{p}_{ij}{c}_j$ which contains the expected value of the residual process at next time conditionally on the occupancy of the state $r_i$ at current time k.

In general, as m varies the model gives predictions for any future time $k+m$.

3. Case Study

The transition towards a sustainable and diversified energy mix is a critical challenge and opportunity for countries worldwide. Italy, with its unique geographical, economic, and social landscape, presents an intriguing case for examining how future energy consumption patterns might evolve. This case study delves into forecasting Italy’s primary energy consumption mix-up to the year 2032, employing the innovative Grey Markov Model MCGM(1,1).

The data for our study were collected from the comprehensive database of "Our World in Data" (https://ourworldindata.org/), renowned for its extensive global energy statistics. This rich dataset provided us with annual figures of primary energy consumption across various sources in Italy, measured in terawatt-hours (TWh), spanning from the year 2000 to 2022. The energy sources detailed in this dataset include Biofuel, Coal, Gas, Hydro, Nuclear, Oil, Solar, and Wind.

Table 1 presents data on primary energy consumption across various sources measured in terawatt-hours for the years 2000 through 2022. Each row represents a year, showing how much energy was consumed from each source during that year. Starting from 2000, when reliance on traditional sources like Coal, Gas, and Oil was predominant, there has been a noticeable diversification and increase in renewable energy sources over the years. For instance, Solar and Wind energy consumption started from nearly negligible amounts in 2000 but show significant growth, indicating a shift towards more sustainable energy sources. Notably, Nuclear energy consumption remains at zero. Following the Chernobyl disaster in 1986, Italy held a referendum in 1987, which resulted in the decision to discontinue the use of nuclear power for energy production. The table is a comprehensive depiction of the evolving energy landscape, highlighting the gradual transition from fossil fuels to renewable energy sources.

Primary energy is the energy that is produced when a source is burned in its raw form. It includes the energy required by the end user in the form of electricity, transportation and heating, as well as inefficiencies and energy lost during the conversion of raw materials into a usable form. Since renewable sources are reported in terms of their electricity output, they are corrected by the ‘substitution method’ for efficiency losses in fossil fuels [29]. It attempts to match non-fossil energy sources to the inputs that would be required if they were generated from fossil fuels. To this end, energy generation from non-fossil sources has been divided by a standard ‘thermal efficiency factor’ by Energy Institute [30]. In this way, it is reliable to compare inefficient fossil fuel inputs to renewable energy sources, which do not have this inefficiency.

In our analysis, we will not use the data referring to biofuel from 2000 to 2003 because they contain zero values which would affect the results unrealistically.

For illustrative purposes, let us revisit all the stages of the previously described model, employing data related to the coal energy consumption in terawatt-hours and set them as the following vector

$$\begin{array}{c}\begin{array}{cc}\hfill X^{\left(0\right)}=(& 145.580,154.697,159.678,172.299,193.013,191.534,193.911,189.848,\hfill \\ \hfill & 183.701,143.823,159.028,178.270,182.733,157.403,153.447,143.395,\hfill \\ \hfill & {127.732,111.904,102.567,77.534,59.253,65.182,84.773)}^T.\hfill \end{array}\end{array}$$

(24)

The vector containing the cumulative values is

$$\begin{array}{c}\begin{array}{cc}\hfill X^{\left(1\right)}=(& 145.580,300.277,459.955,632.254,825.267,1016.801,1210.712,1400.56,\hfill \\ \hfill & 1584.261,1728.084,1887.112,2065.382,2248.115,2405.518,2558.965,2702.36,\hfill \\ \hfill & {2830.092,2941.996,3044.563,3122.097,3181.35,3246.532,3331.305)}^T,\hfill \end{array}\end{array}$$

(25)

and consequently we compute

$$\begin{array}{c}\begin{array}{cc}\hfill Z^{\left(1\right)}=(& 222.9285,380.116,546.1045,728.7605,921.034,1113.7565,1305.636,1492.4105,\hfill \\ \hfill & 1656.1725,1807.598,1976.247,2156.7485,2326.8165,2482.2415,2630.6625,2766.226,\hfill \\ \hfill & {2886.044,2993.2795,3083.33,3151.7235,3213.941,3288.9185)}^T.\hfill \end{array}\end{array}$$

(26)

To estimate the parameters $\widehat{a}$ and $\widehat{b}$, we use the following matrices:

$$\mathit{B}=\left[\begin{array}{cc}-222.9285& 1\\ -380.116& 1\\ ⋮& ⋮\\ -3288.9185& 1\end{array}\right],\mathit{y}=\left[\begin{array}{c}154.697\\ 159.678\\ ⋮\\ 84.773\end{array}\right].$$

(27)

Using Eq. (7) we get the estimation of the developing coefficient $\widehat{a}=0.034$ and the Grey effect $\widehat{b}=210.598$. Meanwhile, we use Eq. (9), to compute the Standard Errors obtaining $SE\left(\widehat{a}\right)=0.006$ and $SE\left(\widehat{b}\right)=13.320$.

Wanting to get a forecast up to 2032, which is ten years ahead of the last available data, we use Eq. (10) to get the cumulated values ${\widehat{x}}_k^{\left(1\right)}$ and Eq. (11) to get the predicted values ${\widehat{x}}_k^{\left(0\right)}$. The results of the GM(1,1) model are shown in Table 2. Moreover, in Figure 1, coal primary and predicted energy consumption in terawatt-hours are plotted. It shows that the prediction starts with the initial available data, rises up, and follows an exponential behaviour.

Now, let us proceed with the correction of predictions through the Markov chain. Consider the initial data and compute the relative errors utilizing Eq. (12). Results are recorded in Table 2. In this series, the minimum is ${\epsilon }_{min}=-23.466$ and the maximum is ${\epsilon }_{max}=80.454$. Now, we split the range of residuals $[-23.466,80.454]$ into $r=4$ intervals and assume each one to be a state

$$\begin{array}{cc}\hfill \mathrm{State}1=& [-23.466,2.514],\hfill \\ \hfill \mathrm{State}2=& (2.514,28.494],\hfill \\ \hfill \mathrm{State}3=& (28.494,54.474],\hfill \\ \hfill \mathrm{State}4=& (54.474,80.454].\hfill \end{array}$$

(28)

In the last column of the Table 2, states are assigned based on the interval each relative error falls into. From Eq. (15), the centers of each state are: $c_1=-11.974$, $c_2=12.490$, $c_3=36.692$, $c_4=69.54$.

Subsequently, the one-step transition probability matrix is computed through Eq. (19)

$$P^{\left(1\right)}=\left[\begin{array}{cccc}0.833& 0.167& 0& 0\\ 0.4& 0.4& 0.2& 0\\ 0& 0.5& 0& 0.5\\ 0& 0.5& 0& 0.5\end{array}\right].$$

(29)

For example, the value $0.833$ represents the probability of reaching state 1 given that the chain at the previous step was in state 1; $0.167$ stands for the probability of reaching state 2 given that the chain was previously in state 1. Raising the matrix to the m-th power makes it possible to obtain the m-steps transition matrix. From the matrix in Eq. (29), we obtain the following stationary distribution

$$\pi =\left[\begin{array}{cccc}0.632& 0.263& 0.053& 0.052\end{array}\right].$$

(30)

Since we know all the real values referred to the years 2000-2022, for these years we use the one-step transition probability matrix updating the state i in Eq. (21). Whereas for the years 2023-2032 we used the state $i=2$, that is the last known state, and the matrix raised to the $m-$th power. Finally, once all the adjustments ${\widehat{\epsilon }}_k^{\left(1\right)}$ are computed, Eq. (20) is used to get the adjusted GM(1,1) predictions from 2001 to 2032. See Table 2 for the results referring to the years of the sample and Table 3 for the prediction from 2023 to 2032. Figure 2 summarizes the coal primary and predicted energy consumption in terawatt-hours through the GM(1,1) model before and after the Marvokian adjustment. The dashed vertical line separates the years of our sample from the years related to the forecast. To its left, as we have previously explained, the one-step transition matrix was used, updating the forecast with the observed value each time, while on the right, the forecast was made starting from the last observed value, using the m-step transition matrix. Compared to the forecast made with the Grey Model, the Markovian correction allows for predictions that follow the trend of the observed data, enhancing the predictive capability.

All the previously described steps were repeated for each energy source, obtaining the predictions from 2023 to 2032. It should be noted that in the cases of Biofuel and Solar, it was not possible to calculate the representative values for all states because some states were never reached, as the relative errors pertained only to certain states. To overcome this problem, we employed Change Point analysis to identify the optimal point at which to truncate each data series [31]. The truncation point was determined to be in the year 2010 for Biofuel and in 2011 for Solar. After implementing these modifications, we were able to determine the representative values for each state and retrace all prior steps. In Table 4 all the GM(1,1) adjusted predictions are shown for all the energy sources while in Table 5 the estimates together with the Standard Errors (provided in parenthesis) of the developing coefficient and the Grey effect for each energy source are provided.

Table 5 contains the data relating to the estimate and the Standard Error of the developing coefficient and the Grey effect for each energy source.

The predictions highlight a steady increase in solar and wind energy, reflecting a global shift towards more sustainable and renewable energy sources. Biofuel predictions also show a consistent yet slight increase, indicating a modest contribution to Italy’s energy mix. Coal consumption, on the other hand, is forecasted to decline, echoing the broader transition away from fossil fuels. The forecast for gas and oil consumption suggests a decrease, aligning with efforts to reduce carbon emissions and dependency on non-renewable energy sources. Hydroenergy consumption shows minor fluctuations, indicating stability in its contribution to the energy mix.

4. Conclusions

In this comprehensive study, we introduced and applied the Grey Markov Model GM(1,1) integrated with Markov Chain Residual Modification to forecast Italy’s power generation mix-up to the year 2032. This approach allowed us to incorporate both the uncertainty inherent in historical data and its fluctuating nature, providing a more accurate prediction model for Italy’s energy future.

Our findings reveal a significant potential for renewable energy sources to increase their share in Italy’s energy mix, indicating a gradual decrease in reliance on fossil fuels. This transition aligns with global efforts towards more sustainable and environmentally friendly energy production methods. Notably, the model demonstrated high accuracy in predictions, affirming the integration of Grey systems and Markov chains as a robust method for forecasting in the energy sector.

The insights garnered from this study underscore the necessity for Italy’s energy policymakers to further encourage and invest in renewable energy sources. Strategic planning should incorporate the predicted shifts in energy production types, emphasizing sustainability and energy independence. Additionally, the model’s predictions can serve as a guide for adjusting current policies and devising new regulations to support the anticipated changes in the energy sector.

While this study offers valuable forecasts and insights, it also opens several avenues for future research:

• Extending the model to other countries and regions to compare and contrast energy mix transitions on a global scale.
• Incorporating additional variables such as technological advancements, policy changes, and economic factors that could impact the power generation mix.
• Exploring the potential impacts of increased renewable energy adoption on grid stability, energy prices, and environmental outcomes.

In conclusion, the application of Grey Markov Models to forecast Italy’s power generation mix presents a promising tool for energy researchers and policymakers alike. By understanding future trends, Italy can better prepare for a sustainable energy future, reducing carbon emissions and fostering a resilient energy system. Further research in this area will undoubtedly refine these predictions and contribute to the global knowledge base on energy forecasting methodologies.