Insurance Contracts for Hedging Wind Power Uncertainty

: This paper presents an insurance contract that the supplier of wind power may subscribe to with an insurance company in order to immunize his/her revenue against the volatility of wind power and prices. Based on empirical evidence, we found that wind power and electricity prices are correlated. Then, we adopted a joint stochastic process to model both time series, which is based on indexed semi-Markov chains for the wind power generation process and on a general Markovian process for the electricity price process. Using a joint stochastic model allows the insurance company to compute the fair premium that the wind power producer has to pay in order to hedge the risk against inadequate revenues. Recursive type equations are obtained for the prospective mathematical reserves of the insurance contract. The model and the validity of the results are illustrated through a real data application.


Introduction
The use of renewable energies is continuously expanding year-by-year, and this trend is expected to persist in the future, although a decline of capacity has been observed related to lockdown measures due to the COVID-19 outbreak during 2020; see [1].
Wind energy covers a large part of renewable energy production even if this industry suffers the high variability of the wind speed process, which in turn impacts the electricity generation.
To remedy this problem, several solutions have been considered. Some of them are of a technical nature; others are based on intangible instruments such as financial products and specific contracts. Among the technical solutions, we remember the recourse to energy storage systems based on batteries that can contribute to absorbing the excess of production during high wind speed periods and to yield required power during low wind speed moments. The literature on this aspect is very large, and we may refer to [2][3][4][5][6] and to the respective bibliographies. Another technical possibility is associated with hybrid power plants that help wind power production with any dispatchable energy source such as gas turbines, air compressed systems, hydroelectric plants, or different; see, e.g., [7]. These engineering-based approaches need sophisticated technical know-how that may be outside of the wind power producer's possibilities. For this reason, alternative solutions based on intangible instruments have been proposed and revealed to be practicable. Some of them are based on the possibility for the wind farm to purchase call/put options to protect itself against the uncertainty of the wind; see [8,9] on devising a barrier option for wind power. A more recent proposal came from [10] where an insurance contract was advanced as an effective tool to hedge wind speed volatility and in turn wind power. The main idea of the insurance-driven approach is to design an insurance contract that may protect the wind power producer (WPP) from undesirable fluctuation of the income generated by the wind power production. The contract is based on an exchange of cash-flows in terms of premium paid and benefits obtained by the WPP that should respect a form of equilibrium. Premiums are paid periodically and are for simplicity constant in time. Benefits are collected in case the wind power sold at the market price of electricity is unable to generate the desired income stipulated in the contract. Accordingly, the paramount objective of this paper is to compute the level of fair premium and the evolution of the mathematical reserves generated by the cash-flows. This approach requires the specification of an accurate wind power model. In [10], the authors applied a weighted indexed semi-Markov chain model (WISMC), which has been proven to possess a suitable probabilistic structure able to reproduce many stylized facts of wind speed and power dynamics; see [11][12][13][14]. This class of stochastic processes is a generalization of semi-Markov processes, which have shown to be very appropriate in a wide spectrum of domains ranging from reliability to finance; see, e.g., [15][16][17].
Nevertheless, the contribution by [10] leaves open space for further improvements. First, it is well known that electricity price evolves randomly in time; see, e.g., [18,19]. Thus, the unrealistic assumption of the deterministic electricity price adopted in [10] should be relaxed in favor of stochastic process assignment. Second, we also know that the electricity price is correlated with the wind power production process because large wind power productions are associated with decreases in the electricity price as an excess of supply on the market; see, e.g., [7,20]. This empirical evidence suggests the adoption of the multivariate model of the wind power generation and electricity price process. For this reason, we propose a multivariate model of wind the energy production and electricity price process using copula functions. After having specified the model, we describe the actuarial structure of the insurance contract, and we derive a recurrent type equation for the expected value of the accumulated discounted benefit function over a given time interval. Consequently, we determine the fair premium and an equation for the prospective mathematical reserve. The model is applied to real data concerning the electricity price process and wind production in Sardinia from 1 August 2014 to 31 July 2018. We show how to fix the fair premium and the computation of the connected actuarial quantities.
From a mathematical point of view, it is relevant to observe that the considered mathematical methods presented in this paper are generalizations of methods proposed in other classical problems of insurance mathematics like life insurance contracts and disability modeling where multi-state models based on Markov chains (see, e.g., [21][22][23]) or semi-Markov chains (see, e.g., [24][25][26][27][28][29][30][31][32][33]) with rewards represent the usual approach. The rest of the paper is organized as follows. In Section 2, we briefly describe the WISMC model of wind energy and the joint model of electricity price and wind power generation. This fixes the framework in which the actuarial evaluation will be subsequently developed. We then present in Section 3 the actuarial investigation with the description of the insurance contract, the computation of the fair premium, and the evolution of the mathematical reserve. In Section 4, we present a numerical example based on real data on wind power production and electricity prices.

The Stochastic Models
In this section, we describe the main stochastic tools adopted for the description of the considered insurance contract. We first present the WISMC model that is used as a marginal model for wind power production. Successively, we introduce a Markovian model for the electricity price process, and using a copula function, we specify a dependence structure between electricity price and wind power production. In this way, we obtain a multivariate model that is at the base of the valuation of the insurance contract.

The WISMC Model of Wind Power Production
Here, we introduce the discrete-time WISMC model with a finite state space in relation to the wind power production of a wind farm.
Indexed semi-Markov models are based on three sequences of random variables defined on a common complete probability space (Ω, F , P). First, we consider the stochastic process (J n , n ∈ N) with values in E = {1, 2, ..., s}. The set E denotes different levels of wind power production. These levels change value in time due to randomness inherent to the wind power production phenomenon, and (J n , n ∈ N) describes the sequences of wind power transitions from one level to another one in time.
Evidently, the wind power production stays at a given level of production for a random time before changing its value. Accordingly, there is another sequence of random variables denoted by (T n , n ∈ N) that indicates the times when the change of states of (J n , n ∈ N) occurs. The sequence (T n , n ∈ N) is an increasing sequence of random times, and their difference (X n := T n+1 − T n , n ∈ N) is called the sequence of sojourn times. Thus, the random variable X n denotes how long the wind power state will remain at level J n before making a transition to level J n+1 .
A detailed description of the system needs the specification of the so-called index process that we are going to introduce.
Assume disposing of a set of past information about wind power levels and transition times, i.e., then consider a generic bounded function f λ : E × N × N → R called the score function, and define the index process as: A specific choice of the score function f λ may depend on the considered application and will be discussed later in the applied section. The quantity λ is the memory parameter and is used to calibrate the overall memory of the process to the real data. Different values of λ could be used according to the need to make more relevant the information brought by the past observations. The calibration of λ will also be discussed later in the application. From an intuitive point of view, the index process I λ n assigns to a past trajectory of wind power production a total score. This is done by attaching to every past wind power observation a score f λ (J n−1−k , T n , a), which depends on the recorded value of the wind power J n−1−k , on when this observation was registered, and on the present time T n .
These three random processes are assumed to interact, and therefore, we specify a dependence structure between them. Specifically, we assume that: where σ(J h , T h , I λ h ), h ≤ n is the natural filtration of the three-variate process. Relation (2) asserts that the knowledge of the last wind power value (J n = i) and of the value of the index process (I λ n = v) suffices to give the conditional distribution of the couple J n+1 , T n+1 − T n whatever the values of the past variables might be. This means that knowing the couple (i, v), we dispose of the joint distribution of the next wind power value J n+1 and of the time when the system will enter into J n+1 . This model is particularly able to increase the memory of the process incorporating past information about the system's evolution. It should be remarked that whenever q λ ij (v; t) should be constant with respect to v for all possible combinations of i, j, t, then the WISMC model collapses into a standard semi-Markov chain model.
Let H λ i (v; ·) be the sojourn time cumulative distribution function in wind power state i ∈ I with index process value v, that is: It expresses the probability to change the actual wind power i in a time less than or equal to t given the indexed process has value v. The standard notation for the survival function is adopted . Let N(t) = sup{n ∈ N : T n ≤ t} be the number of transitions held up to time t, and let Z λ (t) = J N(t) be the state of the system (wind power) at time t. We refer to (Z(t), t ∈ N) as a weighted-indexed semi-Markov chain and to (J n , n ∈ N) as the indexed embedded Markov chain.
The processes {J n , T n , I n (λ)} provide a description of the status of the system when a transition occurs, i.e., in correspondence to any one of the jump times T n . However, during any time t belonging to a sojourn time (that is for T n < t < T n+1 ), while the state of the wind power remains unchanged (J N(t) = J n ), the value of the index process updates due to the information that the embedded Markov chain has not changed state along the sojourn time or, equivalently, the time elapsed since the last transition has increased. The time elapsed since the last transition is called the backward recurrence time process and is defined by B(t) := t − T N(t) . Thus, to completely characterize the status of the system also along a sojourn time, it is necessary to extend the definition of the index process allowing considering any time t ∈ N as follows: where θ = 1 {t>T N(t) } . It is not difficult to realize that if t = T n , then N(t) = N(T n ) = n, and accordingly: Now, observing that the previous relations imply that J N(t) = J n , by substitution in (4), we get: which coincides with (1) and shows that I λ (t) = I λ n . Let us introduce the marginal one-step distribution of WISMC: In our application, it denotes the conditional probability distribution of wind power production at time t + 1 given the values of the wind power production and the backward and index processes at time t.
This probability can be obtained as follows: where The probabilities p (i,u,v) (j, d) can be obtained from the knowledge of the kernel of the WISMC according to Lemma 2.1 and Theorem 2.1 in [34]. We conclude this short overview on the WISMC model by considering two concepts that we need to introduce. They are taken from [35] and play an important role in the determination of the results of this study.
. . , n} be a sequence of states and corresponding transition times, i.e., . . , 0, . . . , n, i α ∈ E, t α ∈ Z}, then we define the shift operator: Intuitively, the shift operator when it is applied to a trajectory (i, t) n+1 −m renders back another trajectory where the sequence of visited states is the same as in the input trajectory with the difference that transition times are translated to t n+1 time units backward and the number of transitions is set one unit backward.
The following assumption concerning the score function f λ will be needed in the rest of the article: Lemma 1. For a WISMC with score function f λ that satisfies Assumption 1, for fixed arbitrary state j and time t and (i, t) n+1 −m ∈ Θ n+1 −m , we have: The result presented in Lemma 1 focuses on a class of score functions leading the probability (7) to be independent of n. Accordingly, the WISMC inherits a homogeneity property that is particularly useful for the applications of the model. Throughout this article, we are going to consider only homogeneous WISMC.

Joint Model of Electricity Price and Wind Power Production
A wind power producer (WPP) sells energy in the spot market. Accordingly, he/she will accept the market price of electricity as a counterpart for his/her power production. The electricity price changes in time randomly, and different stochastic processes have been advanced as a suitable description of the statistical properties of this variable. A very exhaustive review and perspective were given in [18].
In this paper, we focus on a simple structure for the electricity price model that can be merged, successively, with the WISMC model of wind power production. For this reason, we focus on a stochastic process {π e (t), t ∈ N} that satisfies a general Markovian property. Precisely, we assume: Assumption 2. The energy price dynamics satisfies the Markov property, i.e., A large variety of models satisfies this property, and an example is given by: where (t + 1) is a Gaussian noise. Equation (8) defines a first-order autoregressive model (AR(1)) that exhibits a mean reversion property. Extensions to AR(p) models fall well within our framework, as well as processes that allow the presence of a jump component: where η(t + 1) is a pure jump process. The observation of data on electricity prices and wind power production suggests the presence of a dependence structure between these random processes; see, e.g., [7,20]. In general, it is expected that if a large amount of wind power is introduced into a power system, this would tend to lower the spot price on the related energy market. Motivated by this empirical evidence, we advance another assumption: Assumption 3. The wind power production and the electricity price are dependent processes. Their joint probability distribution satisfies the following property: where σ(Z λ (s), π e (s), s ≤ t) is the sigma algebra generated by the WISMC (which includes information on past values of the backward recurrence time process and of the index process as well). The symbol C is used to denote a copula function, which merges the conditional marginal distribution of electricity price: F π e (p; x) = P[π e (t + 1) ≤ x|π e (t) = p], and the conditional marginal one-step distribution of wind power production: Thus, once the marginal models of electricity price and wind power are settled on, the adoption of a specific copula function leads to a joint model that is at the base of the actuarial problem we are going to discuss in the next section.

The Insurance Problem
In this section, we argue that a WPP may protect himself/herself against the variability of the income generated through its activity using actuarial techniques. Let r(t) := Z λ (t) · π e (t) be the WPP's revenues collected at time t by producing a level Z λ (t) of wind power that is sold on the market at the price π e (t). The process r(t) is complex and is subject to consistent variability in time due to its basic components: the produced wind power Z λ (t), the electricity price process π e (t), and their dependence structure.
Here, we consider the possibility to hedge this risk using an insurance contract. At any time s ∈ N, the WPP collects random revenues {r(s)} s∈N . Assume, he/she wishes to guarantee for himself/herself at any time s at least a revenue equal to K ∈ R + . To this end, he/she may sign an insurance contract such that: • if r(s) < K, he/she gets from the insurer the benefit K − r(s); • if r(s) ≥ K, no money transfer from the insurer to the WPP occurs; • at any time during the validity period of the contract, he/she pays a fixed premium to the insurer equal to U ∈ R.
• money amounts are discounted with fixed discount factor v; accordingly, v s denotes the discount factor for s periods of time.
A traditional problem in actuarial mathematics is how to determine the value of the premium U. The simplest way is to fix U in such a way that premiums and benefits are in equilibrium over a considered time interval. Thus, let B(t 0 , u, T) be the total benefits paid by the insurer during the time interval [t 0 + u + 1, T]: The dependence on three time variables is considered with the intention of considering explicitly the effects of the last transition time (t 0 ), the time u elapsed in the last state until current time (t 0 + u), and the horizon time of evaluation T.
Clearly, B(t 0 , u, T) is a random process that inherits statistical properties from the revenue process r(t). A viable solution is to define the total fair premium as the amount U u,T := U · v · ∑ T−u k=0 v k such that U u,T is equal to the expected value of the benefits, i.e., E[B(t 0 , u, T)].
The fair premium can now be determined as the amount U such that: whereä j:T−u| represents the present value of an annuity due, which is a series of unit payments at the beginning of each period for T − u periods discounted using the interest rate j = 1−v v . Thus, the possibility to determine the fair premium passes from the computation of the total expected benefits E[B(t 0 , u, T)]. To this end, denote by B((i, t) 0 −m , u, p; T) the random variable that has the same distribution as the conditional distribution of the random variable B(t 0 , u, T) given that: Finally, denote by: the conditional expectation with respect to the information available up to present time t 0 + u, which includes past values of the wind power production, corresponding transition times, backward recurrence time process value, and electricity price.

Lemma 2.
Under Assumptions 1-3, the accumulated discounted benefit process satisfies the following stochastic relation: and consequently: Proof. To show that: and: share the same distribution, it is sufficient to remark that these random variables are obtained by the summation of the same addends, which are discounted by the same time length. These random addends have the same conditional distribution even if they are defined on different information sets, which are: for (14) and: for (15). Nonetheless, from Assumption A2, we have that: D(π e (t 1 + 1)|π e (t 1 ) = p) = D(π e (1)|π e (0) = p), having denoted by D(X) the distribution of X. Moreover, from Lemma 1 in [35], we can obtain that, given (J, T) 1 −m = (i, t) 1 −m , the index process assumes the value: and the same value is assumed conditional on (J, T) 0 −m−1 = •(i, t) 1 −m , i.e., I λ 0 = x. This implies that: s 0 = i 1 and k 0 = t 1 − t 1 = 0 being from the definition of the shift operator (1). This shows that also the conditional distributions of Z λ (·) are the same in the two cases. The simultaneously validity of Relations (14)- (16) and (18) proves (12) and (13).
The first-order moment of the accumulated discounted benefits can now be calculated using the relation derived in the next theorem. Theorem 1. Under Assumptions 1-3, the expected value of accumulated discounted benefits satisfies the following recurrent equations, for all (i, t) 0 −m ∈ Θ 0 −m , u ∈ N, T ∈ N, u + t 0 < T, and p ∈ R, where t 0 = 0, t 1 = u + 1, and: Proof. Since t 0 = 0, we have that: Let χ(E) be the indicator variable for event E, then the expectation (20) can be further decomposed into: Formula (21) can be obtained according to the integral representation of expectation as follows: Observe now that the event {Z λ (u + 1) = i 0 , Z λ (u) = i 0 , B(u) = u} implies that B(u + 1) = u + 1, and consequently, the former integral becomes The conditional joint distribution of the random vector {Z λ (u + 1), π e (u + 1)} inside Formula (23) can be derived using Assumption 3 as a function of the copula and the marginal distributions of wind power production and electricity price process. Indeed, for small dx, we have: where: and: Firstly, it is worth observing that this joint distribution is of a mixed type, having a continuous part (F π e (p; ·)) and a discrete one (F (i 0 ,u,α 0 ) (·)). Secondly, the expression (25) is in agreement with a general formula for multivariate mixed distributions presented in [36].
The relation: is therefore valid.
To compute the expectation in Formula (21), we observe that the event {Z λ (u + 1) = i 1 = i 0 , Z λ (u) = i 0 , B(u) = u} implies that B(u + 1) = 0. We use then the integral representation of expectation and analogous computations as we did before to get: It remains to compute Formula (20). It is simple to realize that: The latter expression can be rewritten using the tower property of conditional expectation as follows: Observe that the random variable: (K − π e (s)Z λ (s)) + v s−(u+1) |Z λ (u + 1), B(u + 1), π e (u + 1), (J, is distributed the same as: From Lemma 2, we have that: and from this, we get the following representation for Formula (29): By integration, we obtain: Similarly, one obtains for the second part of (31) that: which completes the proof of the theorem.
Remark 1. Theorem 1 can be seen as an extension of Theorem 1 in [29] and also of those obtained in [10] to the more general framework represented by WISMC models. More precisely, our result includes the possibility to manage random rewards that are not only a deterministic function of the state of the chain, but that are obtained through a transformation of a stochastic process (electricity price process), which is dependent through a copula function on the state process (wind power production) modeled using a WISMC model.
Theorem 1 determines an expression for the expected value of paid benefits during the validity period of the contract. Therefore, if the premium is fixed according to Relation (10), at any time s, we have the generation of a payment (positive or negative) equal to: which guarantees a global balance during the time interval [t 0 + u, T]. Anyway, the accumulation process of premiums and benefits is not in equilibrium at any time L ∈ [t 0 + u, T] where an excess of paid premiums or received benefits may occur. To quantify this payment process, which is important to understand the cash-flows associated with the contract, it is sufficient to replace the time T (maturity of the insurance contract) in Equation (19) by time L and to consider payments according to Relation (32).
To this end, if we denote by R(t 0 , u, L) := ∑ L s=t 0 +u+1 ψ(s)v s−u the accumulated rewards from the inception time t 0 + u up to the time L and by R((i, t) 0 −m , u, p; L) the random variable, which has the same distribution as the conditional distribution for the r.v. R(t 0 , u, L) given that the past information −m , u, p; L)] and obtain a recursive equation for the so-called prospective reserve W((i, t) 0 −m , u, p; L), i.e., The quantity W((i, t) 0 −m , u, p; L) expresses the average level of payments up to time L. A positive value of W((i, t) 0 −m , u, p; L) denotes the fact that up to time L, there was an excess of benefits with respect to premiums. Conversely, a negative value of W((i, t) 0 −m , u, p; L) indicates that premiums were paid in excess as compared to the benefits. Obviously, at the maturity, T of the contract premiums and benefits should be in equilibrium; therefore, W((i, t) 0 −m , u, p; T) = 0. The knowledge of the prospective reserve is very important in insurance mathematics because a positive function W((i, t) 0 −m , u, p; ·) highlights the need to charge the fair premium in such a way to avoid that on average, the insurance company will pay in advance. This prevents the contract from not being self-financed.
As a matter of example, we show in Figure 1 two hypothetical paths of the prospective reserve. Both curves start from zero and at the end of the contract assume a null value. The positive curve denotes a contract not self-financed, needing the introduction of a charge over the fair premium due to expected benefits, which forgo expected premiums. The negative curve indicates a contract that is self-financed due to expected premiums forgoing expected benefits.

Materials and Methods
The data covered a 4 year period (from 1 August 2014 to 31 July 2018) on an hourly basis. The wind speed dataset was retrieved from the modern-era retrospective analysis for research and applications, version 2 (MERRA-2) project (https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/).
The data concerning electricity price are freely available from the website of the Italian company "Gestore Mercati Energetici" at (http://www.mercatoelettrico.org/It/download/DatiStorici.aspx). We obtained the data for the Sardinia zone, covering the period described before on an hourly basis.
The time series of wind speed and prices for the analyzed period are shown in Figure 2. We also plotted the empirical probability density function of both variables ( Figure 3) and compared them with a Weibull distribution for wind speed and a lognormal distribution for prices  We also deduced a dependence between price and wind speed that gave rise to low negative correlation that had a small variation during the day; see Figure 4. As a consequence, this correlation should be carefully managed. To obtain the energy produced by one wind turbine, we transformed the wind speed in wind power by means of the power curve of a hypothetical wind turbine. The turbine is located in Sardinia (Italy) and has the following characteristics: To extend the analysis to a wind farm that consists of multiple turbines, as a first approximation, it is possible to multiply the wind energy produced by the number of turbines.
According to Section 2.1, wind power was modeled by using a WISMC model. The parameters were optimized by using the algorithms described in [10,37]. Following the same works, we chose an exponentially weighted moving average (EWMA) of the past wind speed values to define the index process.
Then, wind power was discretized into nine states with the cited optimization algorithm that reduced the mean squared error between the real power autocorrelation function and the simulated one. According to [20], we selected the AR(2) process for the electricity price time series. Then, we modeled correlation between prices and wind power production by using a Gaussian copula. All parameters, for price time series and the copula, were optimized by using the MLE algorithm implemented in MATLAB software. To compare our model with real data, we generated 10 6 synthetic trajectories of the joint process of wind power and electricity prices of the same time length of the real database by means of Monte Carlo simulations.

Results on the Insurance Problem
In Figure 5, we compare the fair premium as defined in Equation (10) for real data and simulated data as a function of the revenue K. We used the hypothesis that the premium and benefits are settled daily. All values were actualized by using a constant interest rate fixed at 1% (the use of a term structure for the interest rate can be implemented easily). In the plot, the dashed line represents simulated data, while the continuous line represents real data. It is obvious from the figure that the model is able to reproduce perfectly the real behavior. As expected, the fair premium is an increasing function of the revenue K.  In Figure 6 the daily difference between benefits and the periodic premium is calculated as a function of time and revenue K. As one would expect, this difference is a random variable, and the variance increases with the revenue value.

K (Euro)
Time (Day) Figure 6. Differences between benefits and the periodic premium as a function of time and revenue K.
As the last result, in Figure 7, we show the prospective reserve, again as a function of time and revenue. According to the theoretical results, the reserve starts from zero and reaches a zero value at the end of the contract. Like the difference between benefits and premium, the prospective reserves deviate from zeros more for increasing values of the revenue K.

Conclusions
In this paper, we showed how to model jointly wind power and electricity prices that, according to real data, have a small negative correlation, which, nonetheless, cannot be ignored. The advanced joint stochastic process is based on indexed semi-Markov chains for the wind power generation process and on a general Markovian process for the electricity price process, which are merged through copula functions. The joint process is used to compute the fair premium and the prospective mathematical reserve for an insurance contract that a wind power producer can stipulate to hedge himself/herself from the risk of low revenues. The results are a generalization to a dependent framework of results obtained in the case of constant electricity prices, which are independent of the power process by [10]. We presented the computation of the most relevant insurance quantities both from a mathematical and an empirical point of view. The empirical results show a very good agreement between real and synthetic data obtained from Monte Carlo simulations. Nevertheless, this mathematical approach suffers from a few limitations. The first is due to its probabilistic structure, which can only be applied in the presence of a consistent dataset with many observations. For example, the same application could not be executed on data sampled on monthly or even weekly time scales. The second limitation concerns the choice of the score function that defines the index process, which is not yet automatized and only done in a heuristic way. Further developments of the model including the computation of the retrospective mathematical reserve and the adoption of more general schemes for computing the premium that go beyond the application of the fair premium principle seem highly promising. Although such considerations on wind insurance are entirely theoretical at the moment, they may have interesting practical implications and allow insurance companies to enter into a new business line.