An assessment of model risk in pricing wind derivatives

Abstract

Wind derivatives are financial instruments designed to mitigate losses caused by adverse wind conditions. With the rapid growth of wind power capacity due to efforts to reduce carbon emissions, the demand for wind derivatives to manage uncertainty in wind power production is expected to increase. However, existing wind derivative literature often assumes normally distributed wind speed, despite the presence of skewness and leptokurtosis in historical wind speed data. This paper investigates how the misspecification of wind speed models affects wind derivative prices and proposes the use of the generalized hyperbolic distribution to account for non-normality. The study develops risk-neutral approaches for pricing wind derivatives using the conditional Esscher transform, which can accommodate stochastic processes with any distribution, provided the moment-generating function exists. The analysis demonstrates that model risk varies depending on the choice of the underlying index and the derivative’s payoff structure. Therefore, caution should be exercised when choosing wind speed models. Essentially, model risk cannot be ignored in pricing wind speed derivatives.

1. Introduction

As the net zero commitments gather global momentum, the wind power industry has experienced record growth in recent years. With almost 94 GW of capacity newly added in 2021, the total global wind power capacity has reached 837 GW (Global Wind Energy Council, 2022). The average annual electricity demand covered by wind power plants increased to 14% in Europe in 2021 (WindEurope, 2022). As more countries seek to meet their renewable energy targets, wind power is expected to comprise a higher share in the electricity market.

Wind turbines depend on the immediate availability of wind for energy production. Given that storing wind for later use is infeasible, the inherently intermittent and non-storable nature of wind results in significant fluctuations in power generation. This variability creates revenue uncertainty for energy producers (Dowell & Pinson, 2016; Zhang et al., 2014). To mitigate the financial impacts of unfavorable wind conditions, wind derivatives have been developed as a risk management tool (Alexandridis & Zapranis, 2013).

The payoffs of wind derivatives are associated with an index representing wind conditions or wind power production in specific geographical areas and time periods. Alexandridis & Zapranis (2013), Benth & Benth (2009), and Rodríguez et al. (2021) suggest using accumulated wind speed over a time period as the underlying index since wind speed is a critical determinant of wind power production. The European Energy Exchange trades wind power futures based on the average load factor, calculated as the ratio of produced wind power to installed production capacity. Nasdaq’s wind power future is written on the German wind index NAREX WIDE, a synthetic index derived from wind conditions and wind turbine power curves.

Wind derivatives make payouts to energy producers when adverse wind conditions result in low wind power production, mitigating revenue losses. This financial protection allows energy producers to maintain stability in their operations and revenues, thereby making renewable energy projects more appealing to investors and promoting the global shift towards cleaner energy sources (Gatzert & Kosub, 2016). By transferring the risk associated with wind power production to financial markets, wind derivatives contribute to the growth and resilience of the renewable energy sector.

As the demand for renewable energy and effective risk management tools grows, further research in the field of wind derivatives becomes increasingly important. In this paper, we contribute to the wind derivative literature by investigating the pricing of wind derivatives, which is fundamental to the development of the wind derivative market.

We consider futures and options written on wind speed indices. The pricing of such wind derivatives necessitates two key components: a wind speed model and a pricing approach. Misspecification of the wind speed model can result in considerable pricing errors. Several researchers (Benth & Benth, 2009; Alexandridis & Zapranis, 2013; Rodríguez et al., 2021) assume normally distributed wind speeds when pricing wind derivatives. However, wind speed data often exhibit skewness and excess kurtosis (Gracianti et al., 2021), violating the normality assumption. This paper aims to examine how misspecification of the wind speed model affects wind derivative prices.

Our contribution to this paper is fourfold. First, we propose the use of the generalized hyperbolic (GHYP) distribution for capturing the skewness and leptokurtosis observed in wind speed data. To the best of our knowledge, this is the first instance of applying the GHYP distribution in wind speed modeling. Our proposed model not only accounts for non-normality but also incorporates various stylized features of wind speed.

Second, we develop risk-neutral pricing approaches that are suitable for both the proposed parametric wind speed model and the semi-parametric model proposed by Gracianti et al. (2021). By applying the conditional Esscher transform for the change of measure, we obtain wind derivative prices.

Third, we compare wind derivative prices under different distribution assumptions and assess how the extent of model risk varies with the underlying index and payoff structure of the derivative. This comprehensive analysis allows us to better understand the impact of distribution assumptions on wind derivative pricing, thereby highlighting the importance of accurate wind speed modeling.

Fourth, the wind speed modeling and wind derivative pricing, as discussed in this paper, exemplify how actuarial expertise can be applied to address the financial uncertainties associated with renewable energy and climate change mitigation. Actuaries are adept at creating sophisticated models to capture the complex dynamics and uncertainties in financial and environmental data. In this paper, actuarial expertise is leveraged to develop more accurate models for wind speed, which is essential for pricing and risk management of wind energy. Actuaries have a strong background in pricing and valuation of complex financial instruments. We demonstrate how actuaries can apply their knowledge to the development of innovative pricing methodologies for wind derivatives, which account for the non-normal characteristics of wind speed data and other factors influencing wind power production. Actuarial science has long been instrumental in quantifying and managing risks, and as the world faces unprecedented challenges due to climate change, the need for actuarial expertise in tackling these risks becomes increasingly urgent. Our work aims to inspire further research in the area of climate-related problems and encourage the development of innovative solutions.

Gracianti et al. (2021) find that wind speed exhibits several stylized features, including temporal correlation, seasonality in both mean and variance, volatility clustering, and non-normality. Although various time series models (Benth & Benth, 2009; Alexandridis & Zapranis, 2013; Rodríguez et al., 2021; Benth & Benth, 2010; Benth & Šaltytė, 2011) have been proposed previously for wind speed, only a subset of the stylized features is considered in these models. Gracianti et al. (2021) capture all these features with a semi-parametric model, in which bootstrap is used to retain the non-normality, and a seasonal-autoregressive-seasonal-generalized autoregressive conditional heteroskedasticity (s-AR-s-GARCH) structure is used to incorporate all other features. The parameters in the s-AR-s-GARCH structure are estimated with the quasi-maximum likelihood (QML) method, which assumes normality. However, the bias of the QML estimator is not studied for the wind speed data.

To reduce potential bias in model estimates, we propose that the GHYP distribution and its two special cases – normal inverse Gaussian (NIG) and hyperbolic (HYP) – be used for capturing the non-normality in wind speed. The GHYP distribution is a flexible distribution class that nests HYP, NIG, Student’s $t$ , variance-gamma, and normal distributions. It has been shown to provide a good fit to daily stock returns, which also exhibit leptokurtosis, by Bibby & Sørensen (2003) and Chen et al. (2008). Our proposed wind speed model amalgamates the s-AR-s-GARCH structure and an error term that follows a GHYP distribution. Although the NIG distribution has been adopted by Ahčan (2012) and Cabrera et al. (2013) for temperature and rainfall modeling, we have yet to see any application of the GHYP distribution in wind speed modeling. The proposed parametric model can be easily estimated using the maximum likelihood (ML) method. By comparing the estimates obtained from the QML method and the ML method with the GHYP distribution, we show that the bias of the QML estimator is significant.

Risk-neutral pricing is widely used for equity derivatives. The idea is to change the real-world measure to a risk-neutral one under which all assets have the same expected rate of return, namely the risk-free rate. There have been multiple attempts to price wind derivatives with the risk-neutral approach. Benth & Benth (2009) represent wind speed dynamics with a continuous-time autoregressive model, and the Girsanov theorem is used to obtain the risk-neutral wind speed dynamics. Alexandridis & Zapranis (2013) employ a similar approach with Benth & Benth (2009) for pricing wind speed futures. Benth & Pircalabu (2018) adopt non-Gaussian Ornstein-Uhlenbeck-type processes for wind power production index and derive the risk-neutral prices of wind power futures using Esscher transform. However, the Girsanov theorem and Esscher transform do not apply to the wind speed dynamics proposed in this paper. The Girsanov theorem is applied to wind dynamics driven by Brownian motions by Benth & Benth (2009). The s-AR-s-GARCH structure with non-normal error terms is unsuitable for the Girsanov theorem. Although Benth & Pircalabu (2018) use a gamma distribution for log wind power index to allow for skewness and leptokurtosis, their model assumes deterministic volatility. The Esscher transform cannot be applied to models with stochastic volatility. Therefore, our proposed wind speed model requires a different approach for the change of measure.

We adopt the conditional Esscher transform to meet our pricing needs. The conditional Esscher transform is developed by Bühlmann et al. (1996) as a generalization of the Esscher transform to serve a more general class of stochastic processes. The conditional Esscher transform has been employed to price derivatives when the underlying asset follows various stochastic processes, including a GARCH process with Gaussian and non-Gaussian innovation distributions (Siu et al., 2004; Badescu & Kulperger, 2008), a regime-switching geometric Brownian motion (Elliott et al., 2005), and a Markov-modulated jump-diffusion model (Elliott et al., 2007). The conditional Esscher transform is well-suited for the proposed wind speed dynamics.

We show two different approaches to using the conditional Esscher transform. In the first approach, the risk-neutral wind speed dynamics are derived using the conditional Esscher transform. The wind derivative prices are determined by simulating wind speed under the risk-neutral measure. In the second approach, wind derivative prices are determined by simulating wind speed under the real-world measure and by evaluating the Radon-Nikodym derivatives associated with the conditional Esscher transform. This approach is also considered by Badescu & Kulperger (2008) for stock option pricing. While the first approach can be more computationally efficient, the second approach is more versatile since it does not require risk-neutral dynamics. For the semi-parametric wind speed model proposed by Gracianti et al. (2021), only the second approach can be applied.

To assess the impact of model misspecification on wind derivative prices, we compare prices obtained under three distribution assumptions: normal, GHYP, and non-parametric. Futures and put options written on two types of wind speed indices are considered in the pricing. We examine whether the extent of model risk changes with the choice of the underlying index and derivative payoff structure.

The remainder of this paper is organized as follows: section 2 describes the fundamentals of wind power generation; section 3 constructs wind speed derivatives; section 4 introduces the wind speed data used in this paper; section 5 presents wind speed models with different distribution assumptions and their estimation; section 6 discusses the risk-neutral pricing of wind speed derivatives using the conditional Esscher transform; section 6 summarizes the pricing results and assesses model risk in wind derivative pricing; and section 7 concludes the paper.

2. Wind Power Generation

In this section, we provide some fundamental facts underlying the operation of wind turbines that are relevant to our research. Modern wind turbines convert the kinetic energy of wind into mechanical energy and subsequently into electrical energy (Manwell et al., 2009). Wind turbines consist of several components, including a rotor with blades, a nacelle housing the generator, and a tower supporting the structure. When wind flows across the turbine blades, it creates lift forces that cause the rotor to rotate, driving the generator to produce electricity.

Wind power generation is influenced by various factors, including wind speed, wind direction, air density, turbine specifications, and geographical location (Manwell et al., 2009). Wind speed is the most crucial factor affecting wind power generation, as the amount of kinetic energy available in the wind is directly related to wind speed. The relationship between wind speed and power output can be described using the power curve of a wind turbine. The power curve, which depends on turbine designs, represents the electrical power generated by a turbine as a function of the wind speed at a specific hub height. There are three important wind speed thresholds in the power curve:

• • Cut-in wind speed, typically around 3–4 m/s: This is the minimum wind speed at which a wind turbine starts generating power.

• • Rated wind speed, typically around 11–16 m/s: The turbine reaches its maximum or rated power output at this wind speed.

• • Cutout wind speed, usually around 20–25 m/s: This is the wind speed at which the turbine shuts down to prevent damage caused by extremely high winds.

Wind turbines start to rotate and generate power when wind speed exceeds the cut-in speed. Power production increases with wind speed until wind speed reaches the rated speed at which the turbine produces the maximum power output. As wind speed continues to increase, the power output remains constant or is reduced slowly. When wind speed exceeds the cutout speed, wind turbines shut down to prevent damage.

3. Wind Speed Derivatives

The payoffs of a wind speed derivative can be associated with an index that measures wind speed or wind power production in pre-specified geographical areas and time periods. For simplicity and transparency, this paper considers two cumulative wind speed indices proposed by Alexandridis & Zapranis (2013) and Caporin & Preś (2012), respectively. The first index is the sum of the daily average wind speed over a pre-specified period, while the second index is the sum of daily average wind speed between the cut-in and cutout speeds over a pre-specified period. Mathematically, these two cumulative wind speed indices (CWSI) over the period of [ $t_0,t_1$ ] are expressed as follows:

(1) \begin{align} \text{CWSI}^{(1)}_{[t_0,t_1]} &= \sum _{t=t_0}^{t_1} W_t, \end{align}

(2) \begin{align} \text{CWSI}^{(2)}_{[t_0,t_1]} &= \sum _{t=t_0}^{t_1} W_t \cdot \unicode{x1D7D9}(l \leq W_t \leq u), \end{align}

where $W_t$ is the average wind speed at the turbine height during period $t$ (we assume daily periods in this paper), $l$ and $u$ are the cut-in and cutout speeds, and $\unicode{x1D7D9}(\cdot )$ is an indicator function which takes the value of 1 if $l \leq W_t \leq u$ and 0 otherwise. Since turbines only operate when wind speed is between the cut-in and cutout speeds, $\text{CWSI}^{(2)}_{[t_0,t_1]}$ is a better indicator of the amount of wind power production. In this paper, we set the cut-in and cutout speeds at 3 m/s and 25 m/s, respectively, which are typical values for wind turbines.¹

We examine wind derivatives with two types of payoff structures: futures and options. Let us denote $D$ as the tick size of the future or option, $K$ as the option strike price, and $F$ as the future price. At the time $t_0$ , the producer enters a derivative contract written on a wind speed index $\text{CWSI}_{[t_0,t_1]}$ as defined in (1) or (2). At the derivative expiration date $t_1$ , the derivative payoff $g(\text{CWSI}_{[t_0,t_1]})$ is expressed as follows:

\begin{align*} \text{Future pay-off:}\,\, &g(\text{CWSI}_{[t_0,t_1]})=D\left (\text{CWSI}_{[t_0,t_1]}-F\right ),\\ \text{Put option pay-off:}\,\, &g(\text{CWSI}_{[t_0,t_1]})=\text{max}\left (D(K-\text{CWSI}_{[t_0,t_1]}),0\right ). \end{align*}

Since a lower wind speed index generally results in lower wind power production, a wind energy producer should short the future or long the put to hedge the uncertainty in wind power production. When the realized wind speed index falls below the future price or strike price, the producer is expected to produce less wind energy and thus suffer from revenue loss. At the same time, the producer receives a positive payoff from the derivatives to offset the revenue loss from wind energy production.

4. The Wind Speed Data

4.1 The historical data

We obtained the hourly average wind speed at the Bergen station in the period of 2017–2021 from the German Weather Service.² The Bergen station is located in Niedersachsen, the state in Germany with the highest number of turbines and installed wind power capacity. We compute the daily average wind speed by taking the mean of the hourly average wind speed on the same day. Missing values are present in the hourly wind speed dataset. If the hourly values are missing for the entire day, the daily average is imputed by the mean of the daily average values on the same day but in other years with data available.

Figure 1. Left panel – time series plot of daily average wind speed at Bergen station in 2017–2021. Right panel – day-of-year average wind speed, with a red LOWESS smoother line.

The left panel of Fig. 1 displays the time series plot of daily average wind speed at the Bergen station, with seasonal patterns clearly visible. The right panel shows the day-of-year average wind speed, calculated as the average wind speed values for the same day over the years 2017–2021, and the red line represents the corresponding LOWESS smoother. The right panel reveals higher daily average wind speed and variability in winter and lower values in summer. Table 1 provides descriptive statistics for the daily average wind speed. The notable skewness (1.39) and excess kurtosis (2.58) indicate a positively skewed and strongly non-normal wind speed distribution.

4.2 Wind speeds at turbine height

The wind speed data at the Bergen station are collected by weather transmitters located 10 m above the ground. Wind turbines, the height of which ranges from under 60 m to over 170 m, experience much stronger wind than weather transmitters. To estimate the wind speed at turbine height, we transform the wind speed obtained from weather stations using the wind profile power law, which is expressed as follows:

(3) \begin{align} W &= V \Big ( \frac{h}{s} \Big )^a, \end{align}

where $W$ is the wind speed at the turbine height $h$ , $V$ is the observed wind speed at the weather station height $s$ , and $a$ is the power law exponent. The power law exponent varies with atmospheric conditions and land features (Spera & Richards, 1979). It is set at 1/7 in several studies (Burton et al., 2011). Burton et al. (2011) suggest the following approximation for the power law exponent:

\begin{align*} a &= \frac{1}{\ln (h/z_0)}, \end{align*}

where $z_0$ is the roughness length representing the surface roughness. Typical roughness length for open farmland with few trees and buildings is $0.03$ m, while typical roughness length for an offshore wind farm is $0.001$ m. In the following analysis, we assume that $s=10$ , $h=90$ , and $z_0=0.03$ , which yield a power law exponent of $a=0.1249$ .

5. The Wind Speed Model

5.1 The s-AR-s-GARCH model

The pricing of wind speed derivatives requires a model that can adequately capture the features of wind speed and provide good wind speed predictions in the contract period. We capture the features of wind speed with the s-AR-s-GARCH model proposed by Gracianti et al. (2021). This model is expressed as follows:

(4) \begin{align} W_t&=\mu _t+y_t \end{align}

(5) \begin{align} y_t &= \sum _{p=1}^{P} a_p y_{t-p}+ \epsilon _t, \end{align}

(6) \begin{align} \epsilon _t &= \sigma _t z_t, \end{align}

(7) \begin{align} z_t & \sim i.i.d.(0,1), \end{align}

where $\mu _t$ is the unconditional mean of the wind speed on day $t$ , $y_t$ is the demeaned wind speed, $\epsilon _t$ is the error term, $P$ is the number of AR terms, $a_p$ is the AR coefficient, $\sigma _t$ is the volatility on day $t$ , and $z_t$ is a random variable following a distribution with zero mean and unit variance.

Table 1. Statistics of the daily average wind speed at the Bergen station in 2017–2021

Seasonality in daily average wind speed and variance of daily average wind speed is incorporated into the model by including the Fourier series in the dynamics of $\mu _t$ and $\sigma _t$ . Volatility clustering is captured by a GARCH $(U,V)$ model for $\sigma _t$ . Mathematically, $\mu _t$ and $\sigma _t$ are further modeled as follows:

(8) \begin{align} \mu _t &= c_0 + \sum _{r=1}^{R} \phi _{c,r} \cos \left (2\pi r\frac{d(t)}{365}\right )+ \sum _{r=1}^{R} \phi _{s,r} \sin \left (2\pi r\frac{d(t)}{365}\right ), \end{align}

(9) \begin{align} \nonumber \sigma _t^2 &= \; \omega +\sum _{r=1}^{S} \gamma _{c,r} \cos \left (2\pi r\frac{d(t)}{365}\right )+ \sum _{r=1}^{S} \gamma _{s,r} \sin \left (2\pi r\frac{d(t)}{365}\right ) \\ &\quad \,+\sum _{u=1}^U\alpha _u \epsilon _{t-u} ^2 + \sum _{v=1}^V\beta _v \sigma _{t-v}^2, \end{align}

where $R$ and $S$ represent the numbers of terms summed in the Fourier series, and $d(t)$ denotes the day of the year and cycles through $1,2,\ldots,365$ .

Gracianti et al. (2021) provide an in-depth discussion of the motivation of this model based on the study of the daily wind speed at various German weather stations. Similar models have also been considered in Taylor et al. (2009). To avoid repetition, we summarize the key arguments for this model presented in Gracianti et al. (2021). As observed in Fig. 1, seasonality presents in both the mean and variance of wind speed. Therefore, the Fourier series is used in both mean and volatility models to capture seasonality. The AR structure is used to incorporate cyclical effects, and the suitability of an AR structure can be verified by examining the ACF and PACF plots of deseasonalized wind speeds. Volatility clustering is found to be significant upon checking the autocorrelation of squared residuals obtained from fitting an AR model to deseasonalized wind speeds. Overall, the s-AR-s-GARCH model demonstrates reasonable goodness of fit and out-of-sample forecasts.

The values of $P$ , $R$ , $S$ , $U$ , and $V$ are selected to minimize the Akaike information criterion (AIC). The considered value ranges are as follows: $P=1,\ldots,5$ , $R=1,2$ , $S=1,2$ , $U=1,2$ , and $V=1,2$ . Gracianti et al. (2021) examined the ACF and PACF plots for the wind speeds at various stations in Germany and found that a maximum of five AR lags is sufficient for most stations. Fourier series with one or two terms have been commonly employed to model seasonality in temperature and wind speed due to its parsimony (Campbell & Diebold, 2005; Taylor et al., 2009; Espen & Jurate, 2012). They have been shown to be sufficient in removing seasonal fluctuations in their wind speed observations. Consequently, we allow the number of Fourier terms, $R$ and $S$ , to take values of 1 or 2. While the GARCH(1,1) model is often sufficient for capturing the volatility dynamics in wind speed data (Taylor et al., 2009; Gracianti et al., 2021), some cases call for a more elaborate GARCH(2,2) structure (Tol, 1997). As a result, ARCH and GARCH orders, $U$ and $V$ , are assumed to take a value of 1 or 2. We further discuss the model estimation and order selection in section 5.2.

5.2 Quasi-maximum likelihood estimation and order selection

Campbell & Diebold (2005) estimate a similar time series model for daily average temperature using the quasi-maximum likelihood (QML) estimation proposed by Bollerslev & Wooldridge (1992). The QML estimation assumes normal distribution and thus maximizes normal log-likelihood. Bollerslev & Wooldridge (1992) motivate the use of the QML estimation by showing that the estimators are not significantly biased when the true data follow a student’s $t$ distribution. Order selection and parameter estimation can be performed simultaneously. We use QML to estimate s-AR-s-GARCH models with different values of $P$ , $R$ , $S$ , $U$ , and $V$ , and we then select the model with the lowest AIC. The AIC is determined as follows:

\begin{align*} \begin{split} AIC &= -2 \ln (L) + 2N_p, \\ \end{split} \end{align*}

where $L$ is the maximized likelihood of the estimated model and $N_p$ is the number of parameters. Based on the maximized quasi-likelihood, the model with $P=2$ and $R=S=U=V=1$ is found to yield the lowest AIC.

The estimated parameters and corresponding standard errors for the selected model are presented in Table 2. All the estimated parameters are significant at 5% significance level, thereby verifying the presence of the identified wind speed features.

Table 2. Estimated parameters and their standard errors for the s-AR-s-GARCH model using QML estimation

5.3 The generalized hyperbolic distribution

Bollerslev & Wooldridge (1992) suggest that the bias and variability of the QML estimator tend to increase with the degree of heteroskedasticity and leptokurtosis. Wind speed also exhibits skewness, the impact of which on the bias of the QML estimator is unclear. Therefore, it is necessary to consider both QML estimation and maximum likelihood (ML) estimation with other distributional assumptions and to investigate how the choice of distribution assumption affects parameter estimates. In addition, several researchers (Benth & Benth, 2009; Alexandridis & Zapranis, 2013) assume normally distributed wind speed when pricing wind derivatives. The alternative distributional assumptions will aid us in assessing pricing errors that may arise from model misspecification.

To capture the skewness and leptokurtosis present in wind speed, we consider the generalized hyperbolic (GHYP) distribution and its two special cases – hyperbolic (HYP) distribution and normal inverse Gaussian (NIG) distribution for $z_t$ . The GHYP distribution is capable of modeling the leptokurtosis exhibited in daily stock returns as shown by Bibby & Sørensen (2003) and Chen et al. (2008). It is a flexible distribution class that nests HYP, NIG, Student’s $t$ , variance-gamma, and normal distributions.

The probability density function (PDF) of a GHYP distribution can be expressed as follows:

(10) \begin{equation} f_{\text{GHYP}}(x;\; \lambda,\alpha,\beta,\delta,\mu ) = \frac{\gamma ^\lambda }{\sqrt{2\pi }\delta ^\lambda K_{\lambda }(\delta \gamma )} \cdot \frac{K_{\lambda -\frac{1}{2}}(\alpha \sqrt{\delta ^2 + (x-\mu )^2})}{(\sqrt{\delta ^2+(x-\mu )^2}/\alpha )^{\frac{1}{2}-\lambda }} \cdot e^{\beta (x-\mu )}, \;\; \, x \in \textbf{R}, \end{equation}

where $\alpha$ , $\beta$ , $\delta$ , and $\mu$ are the shape, skewness, scale, and location parameters, respectively, $\gamma ^2 = \alpha ^2 - \beta ^2$ , and $K_{\lambda }$ (.) is a modified Bessel function. The domain of the parameters are $\mu \in \textbf{R}$ and

\begin{align*} &\delta \geq 0, |\beta |\lt \alpha, \,\text{if}\, \gamma \gt 0;\\ &\delta \gt 0, |\beta |\lt \alpha, \,\text{if}\, \gamma =0;\\ &\delta \gt 0, |\beta |\leq \alpha, \,\text{if}\, \gamma \lt 0. \end{align*}

The moment-generating function of the GHYP distribution can be written as follows:

(11) \begin{equation} M_{\text{GHYP}}(u) = e^{\mu u}\cdot \left ( \frac{\alpha ^2 - \beta ^2}{\alpha ^2 - (\beta + u)^2}\right )^{\frac{\lambda }{2}}\cdot \frac{K_{\lambda }\left (\delta \sqrt{\alpha ^2 - (\beta + u)^2}\right )}{K_{\lambda }\left (\delta \gamma \right )}, \;\; |\beta + u| \lt \alpha. \end{equation}

The expectation and variance of the GHYP distribution are shown below:

\begin{align*} E(X)&=\mu + \frac{\delta \beta }{\gamma } R_{\lambda }(\delta \gamma ),\\ \text{Var}(X)&=\frac{\delta }{\gamma }R_{\lambda }(\delta \gamma ) + \frac{\beta ^2 \delta ^2}{\gamma ^2}S_{\lambda }(\delta \gamma ), \end{align*}

where $R_{\lambda }(.)$ and $S_{\lambda }(.)$ are functions defined as follows:

\begin{equation*} R_{\lambda }(u) = \frac {K_{\lambda +1}(u)}{K_{\lambda }(u)}\,\,\mbox {and}\,\, S_{\lambda }(u) = \frac {K_{\lambda +2}(u)K_{\lambda }(u)-K_{\lambda +1}^2(u)}{K_{\lambda }^2(u)}, \, \textrm{for}\, u \in \textbf {R}^+. \end{equation*}

Since $z_t$ follows a distribution with zero mean and unit variance, we impose the following constraints on $\alpha$ , $\beta$ , $\delta$ , and $\mu$ to standardize the GHYP distribution:

\begin{align*} \mu + \frac{\delta \beta }{\gamma } R_{\lambda }(\delta \gamma ) &= 0,\\ \frac{\delta }{\gamma }R_{\lambda }(\delta \gamma ) + \frac{\beta ^2 \delta ^2}{\gamma ^2}S_{\lambda }(\delta \gamma ) &=1. \end{align*}

Therefore, only three parameters need to be estimated for the GHYP distribution.

The NIG distribution is a special case of the GHYP distribution with $\lambda =-\frac{1}{2}$ . Setting $\lambda$ to $-\frac{1}{2}$ in (10) yields the PDF of the NIG distribution:

\begin{equation*} f_{\text {NIG}}(x;\;\alpha,\beta,\delta,\mu ) = \frac {\gamma ^{-1/2}}{\sqrt {2\pi }\delta ^{-1/2} K_{{-1/2}}(\delta \gamma )} \cdot \frac {K_{-1}(\alpha \sqrt {\delta ^2 + (x-\mu )^2})}{(\sqrt {\delta ^2+(x-\mu )^2}/\alpha )} \cdot e^{\beta (x-\mu )}, \;\; \mbox {for}\, x \in \textbf{R}. \end{equation*}

The HYP distribution is a special case of the GHYP distribution with $\lambda =1$ . The PDF of the HYP distribution is given by

\begin{equation*} f_{\text {HYP}}(x;\;\alpha,\beta,\delta,\mu ) = \frac {\gamma }{\sqrt {2\pi }\delta K_{1}(\delta \gamma )} \cdot \frac {K_{\frac {1}{2}}(\alpha \sqrt {\delta ^2 + (x-\mu )^2})}{(\sqrt {\delta ^2+(x-\mu )^2}/\alpha )^{-\frac {1}{2}}} \cdot e^{\beta (x-\mu )}, \;\; \mbox {for} \, x \in \textbf {R}. \end{equation*}

5.4 Maximum likelihood estimation

We employ the maximum likelihood (ML) method to estimate s-AR-s-GARCH models with the NIG, HYP, and GHYP distributions, considering various values of $P$ , $R$ , $S$ , $U$ , and $V$ . Subsequently, we select the model yielding the lowest AIC under each distribution assumption. The parameters $\alpha$ and $\beta$ in the NIG, HYP, and GHYP distributions are not estimated directly. We estimate their transformations, $\rho =\beta/\alpha$ and $\zeta =\delta \sqrt{\alpha ^2-\beta ^2}$ , instead. We note that there is a potential identification problem in fitting the model with the GHYP distribution since a different combination of parameters $\lambda$ , $\zeta$ , and $\rho$ may lead to the same or close likelihood.

For all three distributions, the model with $P=3$ and $R=S=U=V=1$ leads to the lowest AIC. It is important to note that, under the assumption of a normal distribution, the optimal AR order $P$ is 2. This deviation in optimal AR order emphasizes how distribution assumptions can impact the model selection, which in turn affects future wind speed predictions, exemplifying model risk.

Table 3 presents the parameter estimates and their standard errors under four different distribution assumptions. The QML estimates are listed under the column “Normal” for comparative purposes. The models with the NIG, HYP, and GHYP distributions have very close parameter estimates since the NIG and HYP distributions are special cases of the GHYP distribution. All three distributions allow for skewness and excess kurtosis. Due to the different values of $\lambda$ used for the three distributions, the estimates for $\rho$ and $\xi$ change with the assumed distribution.

Table 3. Parameter estimates for the s-AR-s-GARCH model assuming various distributions

Let us compare the parameter estimates obtained under the normal and HYP assumptions. Despite differing AR orders under normal and HYP assumptions, parameter estimates remain comparable due to the utilization of the same model structure. Among the five common parameters – $c_0$ , $a_1$ , $a_2$ , $\phi _{s,1}$ , and $\phi _{c,1}$ – used for the conditional mean, only the AR coefficients $a_1$ and $a_2$ exhibit a significant difference under the two assumptions. In our model, the AR coefficients quantify how past shocks influence the current value, while the constant $c_0$ and the Fourier series coefficients $\phi _{s,1}$ and $\phi _{c,1}$ determine the long-term mean of wind speed. Since the estimated values for $c_0$ , $\phi _{s,1}$ , and $\phi _{c,1}$ do not display a significant difference under the two assumptions, we anticipate that the mean forecast of wind speed over the long term will not vary considerably. However, short-term forecasts may show significant differences.

The estimates for all five parameters – $\omega$ , $\alpha _1$ , $\beta _1$ , $\gamma _{s,1}$ , and $\gamma _{c,1}$ – in the conditional volatility display substantial differences under the two assumptions. Furthermore, the estimates for the shape and skewness parameters under the HYP distribution are significant, suggesting a considerable deviation from a normal distribution. Taking these observations into account, it is expected that the shapes of the distribution of future wind speeds under the two assumptions will differ significantly.

Therefore, the choice of distribution assumption impacts both the selected model and the estimates of model parameters, particularly those for volatility, indicating potentially significant model risk. We should be cautious about using the QML estimates in any study where volatility is an important consideration.

5.5 Model comparison

To compare the goodness of fit of the four estimated models, we compute their AIC values. A lower AIC indicates better goodness of fit. Table 4 summarizes the AIC values for the four estimated models. The normal distribution results in a much higher AIC than any other considered distribution, indicating that it is inadequate for the wind speed data. The HYP distribution leads to the lowest AIC, but it only outperforms the NIG and GHYP distributions marginally.

Table 4. AIC values of the four models with different distribution assumptions

To further examine the suitability of the four distributions, we compare the kernel densities of the standardized residuals obtained from the four models in Fig. 2. All four density curves present a significant positive skewness. The kernel density curves obtained under the NIG, HYP, and GHYP distributions overlap with each other. Also, they are slightly more skewed than the density curve obtained under the normal distribution. The overall difference between the four density curves is small. Therefore, the choice of distribution assumption has little influence on the standardized residuals.

Figure 2 Kernel densities of the standardized residuals under various distribution assumptions.

Under each distribution assumption, we also compare the theoretical or estimated density with the kernel density of the standardized residuals in Fig. 3. The solid line in the top left panel represents the theoretical density of the assumed standard normal distribution. In contrast, the dashed line represents the kernel density of the standardized residuals obtained under the normal assumption. The theoretical density is symmetric, while the kernel density is positively skewed. The kernel density also appears to have a heavier right tail and a lighter left tail than a normal distribution. The large discrepancy between the theoretical and empirical distributions of the standardized residuals indicates that the normal assumption is inappropriate. Under the NIG, HYP, and GHYP distributions, the estimated densities are calculated using the estimated distribution parameters: $\lambda$ (only estimated for the GHYP distribution), $\alpha$ , $\beta$ , $\mu$ , and $\delta$ . The differences between the estimated and kernel densities of the standardized residuals are relatively small, suggesting that the NIG, HYP, and GHYP distributions are appropriate for the wind speed data.

Figure 3 Theoretical/estimated densities and kernel densities of the standardized residuals under various distribution assumptions.

5.6 Simulating wind speed from the parametric models

Let us represent the last day in the historical wind speed data by $t=0$ . Given the choice of distribution and the corresponding estimated s-AR-s-GARCH model, we can simulate a path of future wind speed using the following procedures for $t=1,2,\ldots,t_1$ sequentially:

1. i. Calculate the conditional volatility $\sigma _t$ using (9).

2. ii. Simulate random values for $z_t$ from a standard normal distribution if normality is assumed or from the estimated distribution if one of NIG, HYP, and GHYP distributions is chosen.

3. iii. Determine the simulated wind speed value $W_t$ using (4)–(8).

4. iv. Convert the wind speed to turbine height using Equation (3).

5.7 Simulating wind speed from a semi-parametric model

The main reasons for using QML estimation are that the true distribution of the data is not required and that the normal assumption leads to an easy and efficient estimation. However, the normality assumption does not allow us to incorporate skewness and excess kurtosis in the simulation of future wind speed. Such features, which have been observed in Fig. 2, can potentially impact derivative prices significantly. To overcome this problem, we consider the semi-parametric model for wind speed proposed by Gracianti et al. (2021). A non-parametric distribution is used for $z_t$ in the semi-parametric approach, such that the skewness and leptokurtosis can be retained in the simulated wind speed.

Figure 4 Kernel densities of simulated wind speed at turbine height on 03/31/2022 under various distribution assumptions.

The standardized residuals are first obtained from the QML estimation. Sequential procedures similar to those described in section 5.6 are used for simulating future wind speed, but with the simulated normal random values replaced by the bootstrapped samples from the standardized residuals. The simulated values for volatility and wind speed are then calculated using the QML estimates. More specifically, the simulation procedures are summarized as follows:

1. i. Obtain the standardized residuals from the QML estimation of the s-AR-s-GARCH model.

2. ii. Draw a random sample from the standardized residuals and use it as the simulated values for  ${z}_t$ .

3. iii. Calculate the volatility $\sigma _t$ using Equation (9) and the QML estimates.

4. iv. Determine the simulated wind speed value $W_t$ using Equations (4)–(8) and the QML estimates.

5. v. Convert the wind speed to turbine height using Equation (3).

A similar semi-parametric approach is adopted by Badescu & Kulperger (2008), who approximate the true distribution of stock returns with a non-parametric density estimator. They find that their proposed semi-parametric approach yields relatively accurate option prices. Although the semi-parametric approach retains the non-normality in the simulation of future wind speed, it still suffers from biased estimates due to the use of the QML estimator. We investigate the extent to which the semi-parametric model can alleviate the impact of model risk on wind derivative pricing.

5.8 Simulated future wind speed

Fig. 4 plots the kernel densities of the simulated wind speed on 03/31/2022 at the turbine height using models with different distribution assumptions. The vertical lines represent corresponding mean values. The density curves and mean values obtained under the NIG, HYP, and GHYP assumptions almost overlap in the right panel. This observation is not surprising since the three assumptions yield very similar parameter estimates. In the left panel, the density curves under the bootstrap, normal, and HYP assumptions have notably different spreads, with the largest spread under the normal assumption and the lowest spread under the HYP assumption. While similar tail shapes are observed under the bootstrap and HYP assumptions, the normal assumption leads to a lighter right tail and a heavier left tail than the other assumptions. The mean values under the normal and HYP assumptions are very close, while the mean value under the bootstrap method is slightly higher. These observations align with the findings from the comparison of estimated model parameters. The simulated wind speed for 03/31/2022 represents a 90-day-ahead prediction, thereby indicating a long-term projection. Based on the analysis of model estimates, it is anticipated that the mean forecasts for the long-term will be similar under the two distribution assumptions, while the shapes of the distributions exhibit substantial differences.

The observations from Fig. 4 suggest that if the wind derivative prices are largely determined by the mean values of future wind speed, the pricing error due to the normal assumption can be small and the semi-parametric approach may not provide much benefit. However, if the derivative prices are related to a specific portion of the density curve, such as the tails, the pricing error arising from an inappropriate distribution assumption can be large, and simulating from the semi-parametric model may reduce the error to some extent.

Since the NIG, HYP, and GHYP distributions result in very similar parameter estimates and simulated future wind speed, we only consider the HYP distribution, which leads to the lowest AIC in model fitting, in the following analysis.

6. Risk-Neutral Pricing with the Conditional Esscher Transform

6.1 Risk-neutral pricing and the Esscher transform

The wind speed dynamics described in section 5 are specified under the real-world measure. The derivative price can be calculated as the expected discounted derivative payoff under the real-world measure. Since investors demand risk premiums for bearing the uncertainty, the discount rate needs to be adjusted for individual risk preferences; these are difficult to quantify. Risk-neutral pricing adjusts the probabilities of future outcomes such that all assets have the same expected rate of return, the risk-free rate, under the new probability measure. Such a probability measure is called the risk-neutral measure. The price of a derivative is its expected discounted payoff in the risk-neutral world, with the discount rate set to the risk-free rate.

In a complete market where derivative payoffs can be replicated by a portfolio of existing assets, the risk-neutral measure is unique. However, the wind derivative market is incomplete partly because the wind speed is not directly tradable and partly because the sophisticated dynamics of wind speed imply market incompleteness. The incompleteness of the wind derivative market leads to infinite risk-neutral measures, each of which produces a risk-neutral price.

The Esscher transform is a technique used to obtain a risk-neutral measure. The use of the Esscher transform has been justified by several researchers. Gerber & Shiu (1994) show that Esscher transform maximizes the expected power utility function of a representative economic agent. Chan (1999) demonstrates that the minimum relative entropy measure can be constructed from the Esscher transform. Wang (2003) derives the Esscher transform from the Bühlmann’s equilibrium pricing model under some assumptions on the aggregate risk. Kijima (2006) proves that multivariate Esscher and Wang transforms coincide when the underlying risks are normally distributed. Since the introduction of the Esscher transform by the seminal work of Gerber & Shiu (1994), the Esscher transform has been extensively applied in derivative pricing.

6.2 Conditional Esscher transform

There have been several applications of Esscher transform in pricing weather derivatives, including temperature derivatives (Ahčan, 2012), rainfall futures (Cabrera et al., 2013), and wind power futures (Benth & Pircalabu, 2018). The Esscher transform is well-suited for applications that assume deterministic volatility in weather variables. However, our wind speed model incorporates a GARCH volatility structure, necessitating the use of the more versatile conditional Esscher transform. This was proposed by Bühlmann et al. (1996) as a generalization of the Esscher transform to accommodate a broader range of stochastic processes.

Consider a discrete-time economy with the time index $t = 0, 1,\ldots,T$ . Let $(\Omega,\mathcal{F},\mathcal{F}_t,P)$ be a filtered probability space, where $P$ is the real-world probability measure and the filtration $\mathcal{F}_t$ is a sequence of increasing $\sigma$ -fields of $\mathcal{F}$ representing all the market information up to time $t$ . We assume that $\mathcal{F}_0=\{0,\Omega \}$ and $\mathcal{F}_T=\mathcal{F}$ .

Let $(x_t)_{0\leq t\leq T}$ be a process adapted to the filtration $\mathcal{F}_t$ . The conditional moment-generating function of $x_t$ w.r.t. $\mathcal{F}_{t-1}$ is expressed as follows:

\begin{equation*}M^P_{x_t|\mathcal {F}_{t-1}}(u) = E^P[e^{u x_t}|\mathcal {F}_{t-1}].\end{equation*}

We assume that the $M^P_{x_t|\mathcal{F}_{t-1}}(u)$ exists for $0\leq t \leq T$ and $|u|\lt k$ , where $k$ is some positive real number.

Let the process $Z_t$ be defined by

\begin{equation*} Z_t = \prod _{k=1}^t \frac {e^{\theta _k x_k}}{M^P_{x_k|\mathcal {F}_{k-1}}(\theta _k)}, \end{equation*}

where $Z_0=1$ and $\theta _t$ is a $\mathcal{F}_{t-1}$ -measurable process. $\theta _t$ is also called the conditional Esscher parameter. It is easy to verify that the sequence $\{Z_t\}$ is a martingale. The conditional Esscher transform $Q$ w.r.t. $P$ of the process $x_t$ is defined as follows:

(12) \begin{equation} \frac{\text{d} Q}{\text{d}P}\Big |_{\mathcal{F}_t}= Z_t. \end{equation}

The conditional Esscher transform is very versatile and can be applied in the GARCH framework with any type of distribution as long as the moment-generating function exists.

The conditional moment-generating function of $x_t$ under $Q$ can be written as follows: (Siu et al., 2004):

(13) \begin{equation} M^{Q}_{x_t|\mathcal{F}_{t-1}}(u) = \frac{M^{P}_{x_t|\mathcal{F}_{t-1}}(u+\theta _t)}{M^{P}_{x_t|\mathcal{F}_{t-1}}(\theta _t)}. \end{equation}

Equation (13) will be used to derive the wind speed dynamics under the risk-neutral measure.

6.3 The conditional Esscher parameter

In the option pricing literature (Siu et al., 2004; Badescu & Kulperger, 2008) where $x_t$ is the log stock return, $\theta _t$ is the unique solution to the following equation:

\begin{equation*} M^{P}_{x_t|\mathcal {F}_{t-1}}(1+\theta _t)=e^{r}M^{P}_{x_t|\mathcal {F}_{t-1}}(\theta _t), \end{equation*}

where $r$ is the risk-free interest rate. This equation is equivalent to $E^{Q}[e^{x_t}|\mathcal{F}_{t-1}]=e^r$ , which ensures that the expected return of the stock under $Q$ is equal to the risk-free rate. However, we cannot impose this condition when $x_t$ represents wind speed since wind speed is not a tradable asset.

In this paper, we assume that the conditional Esscher parameters are time-invariant and $\theta _t=\theta$ . We examine the prices of wind speed derivatives assuming various values of $\theta$ . In future research, we plan to calibrate the conditional Esscher parameters using the observed wind power future prices and study how the parameters evolve over time.

The absolute value of $\theta _t$ can also be viewed as the market price of risk. The higher the absolute value, the more compensation investors demand for bearing the risk. The sign of $\theta _t$ is not necessarily positive. Let us consider a wind energy producer who trades wind speed derivatives to hedge the uncertainty in wind power production. To achieve an effective hedge, the producer should take a short position on the underlying wind speed index. Therefore, the producer should short the future or long the put option. Since the investors who are in long positions on the underlying index demand risk premiums, the conditional Esscher parameters should be negative such that the distribution of the wind speed index is shifted to the left with the change of measure from $P$ to $Q$ .

6.4 Wind speed dynamics under the risk-neutral measure

Under the real-world measure $P$ , the wind speed is assumed to follow the dynamics expressed in Equations (4)–(9). The mean and variance of the conditional variable $W_t|\mathcal{F}_{t-1}$ are $\mu _t+\sum _p a_p y_{t-p}$ and $\sigma _t^2$ respectively. Assuming that $z_t$ follows a parametric distribution, Normal, NIG, HYP, or GHYP, we derive the distribution of $W_t|\mathcal{F}_{t-1}$ under $Q$ and the corresponding risk-neutral wind speed dynamics.

Normal distribution

Assume that $z_t$ follows a standard normal distribution under $P$ . We have ${W}_t |\mathcal{F}_{t-1} \sim N(A_t,\sigma _t^2)$ , where $A_t=\mu _t+\sum _p a_p y_{t-p}$ . The moment-generating function of ${W}_t |\mathcal{F}_{t-1}$ under $P$ can be written as follows:

\begin{equation*}M_{{W}_t |\mathcal {F}_{t-1}}^P(u) =e^{A_t u+\frac {1}{2}\sigma _t^2 u^2}.\end{equation*}

Using Equation (13), the moment-generating function of ${W}_t |\mathcal{F}_{t-1}$ under $Q$ can be derived as follows:

\begin{align*} M^{Q}_{W_t|\mathcal{F}_{t-1}}(u) &= \frac{M^{P}_{W_t|\mathcal{F}_{t-1}}(u+\theta _t)}{M^{P}_{W_t|\mathcal{F}_{t-1}}(\theta _t)} \\ &=\frac{e^{A_t (u+\theta _t) + \frac{1}{2} \sigma _t^2 (u+\theta _t)^2}}{e^{A_t \theta _t + \frac{1}{2}\sigma _t^2 \theta _t^2}}\\ &= e^{A_t u + \frac{1}{2} \sigma _t^2 u^2 + \sigma _t^2 \theta _t u}\\ &= e^{(A_t + \sigma _t^2\theta _t )u + \frac{1}{2} \sigma _t^2 u^2}. \end{align*}

Therefore, the distribution of $W_t|\mathcal{F}_{t-1}$ under $Q$ is normal with mean $A_t + \sigma _t^2\theta _t$ and variance $\sigma _t^2$ .

Let $S_t=\sum _{r=1}^{S} \gamma _{c,r} \cos \left (2\pi \frac{d(t)}{365}\right )+ \sum _{r=1}^{S} \gamma _{s,r} \sin \left (2\pi \frac{d(t)}{365}\right )$ . The volatility dynamics under $P$ can be rewritten as follows:

\begin{equation*} \sigma _t^2 = \omega + S_t+ \alpha _1 \epsilon _{t-1} ^2 + \beta _1 \sigma _{t-1}^2. \end{equation*}

It is important to note that the distribution of $\epsilon _t|\mathcal{F}_{t-1}$ under $Q$ is normal with mean $\sigma _t^2 \theta _t$ and variance $\sigma _t^2$ . We define a new random variable $\xi _t=\epsilon _t-\sigma _t^2\theta _t$ . The distribution of $\xi _t|\mathcal{F}_{t-1}$ is normal with mean 0 and variance $\sigma _t^2$ under $Q$ . Replacing $\epsilon _t$ with $\xi _t+\sigma _t^2\theta _t$ , the dynamics of wind speed under $Q$ can be written as follows:

(14) \begin{align} W_t&=A_t+ \sigma _t^2 \theta _t+\xi _t, \end{align}

(15) \begin{align} \sigma _t^2& = \omega +S_t+ \alpha _1 (\xi _t+ \sigma _t^2\theta _t) ^2 + \beta _1 \sigma _{t-1}^2, \end{align}

where $\xi _t|\mathcal{F}_{t-1}\sim N(0,\sigma _t^2)$ .

GHYP distribution

Assume that $z_t$ follows $\text{GHYP}(\lambda,{\alpha },{\beta },{\delta },{\mu })$ under $P$ . The conditional variable $W_t|\mathcal{F}_{t-1}$ also follows a GHYP distribution as shown below:

\begin{equation*}W_t|\mathcal {F}_{t-1}\sim \text {GHYP}(\lambda,\frac {\alpha }{\sigma _t},\frac {\beta }{\sigma _t},{\delta }{\sigma _t},A_t+\mu \sigma _t).\end{equation*}

Using Equation (11), the moment-generating function of $W_t|\mathcal{F}_{t-1}$ under $P$ can be written as follows:

\begin{equation*}M^P_{W_t|\mathcal {F}_{t-1}}(u) = e^{(A_t+{\mu }\sigma _t) u} \cdot \Bigg ( \frac {{\alpha }^2 - {\beta }^2}{{\alpha }^2 - ({\beta } + u\sigma _t)^2}\Bigg )^{\frac {\lambda }{2}}\cdot \frac {K_{\lambda }({\delta } \sqrt {{\alpha }^2 - ({\beta } + u\sigma _t)^2})}{K_{\lambda }({\delta } \sqrt {{\alpha }^2 - {\beta }^2}) }, \;\; |{\beta } + u\sigma _t| \lt {\alpha }. \end{equation*}

Using Equation (13), the conditional moment-generating function under $Q$ can be derived as follows:

\begin{align*} \begin{split} M^{Q}_{W_t|\mathcal{F}_{t-1}}(u) &= \frac{M^{P}_{W_t|\mathcal{F}_{t-1}}(u+\theta _t)}{ M^{P}_{W_t|\mathcal{F}_{t-1}}(\theta _t)} \\ &=e^{(A_t+\mu \sigma _t)u}\cdot \Bigg ( \frac{{\alpha }^2 - ({\beta } + \theta _t\sigma _t)^2}{{\alpha }^2 - ({\beta } + (u+\theta _t)\sigma _t)^2}\Bigg )^{\frac{\lambda }{2}} \cdot \frac{K_{\lambda }({\delta } \sqrt{{\alpha }^2 - ({\beta }+(u+\theta _t)\sigma _t)^2})}{K_{\lambda }({\delta } \sqrt{{\alpha }^2 - ({\beta } + \theta _t\sigma _t)^2})}. \end{split} \end{align*}

Therefore, under $Q$ , the conditional variable follows a GHYP distribution as shown below:

\begin{equation*} W_t|\mathcal {F}_{t-1}\sim \text {GHYP}\!\left(\lambda,\frac {\alpha }{\sigma _t},\frac {{\beta }+\theta _t \sigma _t}{\sigma _t},{\delta }\sigma _t,A_t+{\mu }\sigma _t\right).\end{equation*}

We define a new random variable $\xi _t$ , which takes the following conditional distribution:

\begin{equation*}\xi _t|\mathcal {F}_{t-1}\sim \text {GYHP}\!\left(\lambda,\frac {\alpha }{\sigma _t},\frac {{\beta }+\theta _t \sigma _t}{\sigma _t},{\delta }\sigma _t,{\mu }\sigma _t\right)\end{equation*}

under $Q$ . It is easy to see that $\epsilon _t|\mathcal{F}_{t-1}$ follows the same distribution with $\xi _t|\mathcal{F}_{t-1}$ under $Q$ . The dynamics of the wind speed under $Q$ can be written as follows:

(16) \begin{align} W_t&=A_t+ \xi _t, \end{align}

(17) \begin{align} \sigma _t^2& = \omega +S_t+ \alpha _1 \xi _t ^2 + \beta _1 \sigma _{t-1}^2. \end{align}

NIG and HYP distributions

NIG and HYP are special cases of GHYP distribution. The dynamics of the wind speed under $Q$ are the same as those shown in Equations (16) and (17) with the GHYP parameter $\lambda$ set to $-1/2$ for the NIG distribution and 1 for the HYP distribution.

6.5 Pricing approaches

Under the risk-neutral measure, the price of a wind speed derivative is calculated as the expectation of derivative payoffs discounted at the risk-free rate. Therefore, the derivative price at time 0 can be written as $E^Q[e^{-rt_1}g(\text{CWSI}_{[t_0,t_1]})]$ . We consider two approaches for the expectation calculation, with one approach simulating future wind speed paths under the risk-neutral measure and the other approach under the real-world measure.

In the first approach, we simulate $N$ paths of future wind speed directly from the risk-neutral dynamics, which are defined in Equations (14) and (15) under the normal assumption and in Equations (16) and (17) under the NIG/HYP/GHYP assumption. The simulation procedures are the same as those described in section 5.6, but they use the risk-neutral wind speed dynamics. Let $\text{CWSI}^{Q}_{[t_0,t_1]}(i)$ denote the cumulative wind speed index determined under the $i$ th simulated risk-neutral path. The derivative price can be calculated as follows:

(18) \begin{equation} \frac{1}{N}\sum _{i=1}^{N}e^{-rt_1} g\left (\text{CWSI}^Q_{[t_0,t_1]}(i)\right ). \end{equation}

This approach is efficient and convenient to use when the risk-neutral dynamics can be easily simulated. It is commonly used in stock option pricing (Siu et al., 2004).

In the second approach, we simulate $N$ paths of future wind speeds from the real-world dynamics as defined in Equations (4)–(7). Let $\text{CWSI}^{P}_{[t_0,t_1]}(i)$ denote the cumulative wind speed index corresponding to the $i$ th simulated real-world path. The derivative price at time $t$ can be calculated as follows:

(19) \begin{equation} \frac{1}{N} \sum _{i=1}^N e^{-rt_1}g\left (\text{CWSI}^P_{[t_0,t_1]}(i)\right ) \frac{dQ}{dP}(i), \end{equation}

where $\frac{dQ}{dP}(i)$ is the Radon-Nikodym derivative defined in Equation (12) with $t=t_1$ and evaluated under the $i$ th real-world path. This approach does not require the risk-neutral dynamics of $W_t$ . It is particularly useful when the risk-neutral dynamics follow a complex structure and thus cannot be simulated easily.

In section 5.7, we describe the simulation of future wind speed from a semi-parametric model, which aims to mitigate potential pricing errors arising from model misspecification. Due to the use of the bootstrap procedure, we cannot express the corresponding risk-neutral dynamics easily. Therefore, only the second pricing approach can be used with the semi-parametric model. After simulating each wind speed path under the real-world measure, we calculate the Radon-Nikodym derivative under this path. Equation (19) is used to calculate the derivative price after simulating all $N$ paths.

7. Pricing Results and Model Risk Assessment

7.1 Simulated wind speed index values under $Q$

Benth & Pircalabu (2018) calibrate the Esscher parameter $\theta$ by utilizing observed prices of wind power futures based on the NAREX WIDE index from 2016. Their findings reveal that the magnitude of the calibrated Esscher parameters ranges between 0 and 0.1 and decreases as the delivery period lengthens. For a three-month contract, the calibrated $\theta$ is approximately 0.6. Furthermore, Benth & Pircalabu (2018) note that the closing price falls below the theoretical price when $\theta$ equals 0, signifying that the energy producer, who holds the short position, pays a positive risk premium. To account for this observation, the condition $\theta \leq 0$ is imposed in our pricing applications.

We select the three values for $\theta$ : 0, $-0.05$ , and $-0.1$ , to illustrate the pricing of wind speed derivatives. We simulate 10,000 wind speed paths for the period of 01/01/2022–03/31/2022 under the risk-neutral measure with $\theta$ set to 0, $-0.05$ , and $-0.1$ , respectively. We view 12/31/2021 as day 0 and assume that the derivatives are traded at $t=0$ . The wind speed index is accumulated over the period of $[t_0,t_1]$ with $t_0=1$ and $t_1=90$ .

Fig. 5 presents the density curves of the simulated $\text{CWSI}_{[t_0,t_1]}$ under various distribution assumptions and values of $\theta$ . The top panels compare the kernel densities of simulated $\text{CWSI}_{[t_0,t_1]}^{(1)}$ while the bottom panels compare those of simulated $\text{CWSI}_{[t_0,t_1]}^{(2)}$ . When $\theta =0$ , simulating under $Q$ is the same as simulating under $P$ . In the top left and bottom left panels where $\theta =0$ , we compare the density curves of three sets of index values simulated under $P$ using the normal, HYP, and non-parametric (bootstrap) distributions, respectively. The vertical lines represent the mean values corresponding to the densities. As $\theta$ becomes more negative, the density curves and their mean values shift further to the left, thereby resulting in lower future prices and higher put option prices.

Figure 5 The kernel densities of simulated $\text{CWSI}_{[t_0,t_1]}$ under various distribution assumptions and $\theta$ values, with the vertical lines representing corresponding mean values.

The spread of the density curve under the HYP distribution also changes with $\theta$ , while the spread under the normal distribution does not experience visible change. The change of measure in risk-neutral pricing shifts the skewness parameter of the HYP distribution by $\theta$ and subsequently changes the variance of the wind speed. As $\theta$ becomes more negative, the variance of wind speed under the HYP distribution is further decreased. On the other hand, the change of measure only shifts the mean but not the variance of the normal distribution.

The top panels show that the mean values of simulated $\text{CWSI}_{[t_0,t_1]}^{(1)}$ under the normal and HYP distributions are very close, although the variance under the HYP distribution is smaller. As shown in Fig. 4, we observed similar mean values of simulated wind speed under the normal and HYP distributions in the real world. Since $\text{CWSI}_{[t_0,t_1]}^{(1)}$ is simply the sum of wind speed values over the period of $[t_0,t_1]$ , the mean values of simulated $\text{CWSI}_{[t_0,t_1]}^{(1)}$ under the two distribution assumptions are expected to be close too. The density curve using the bootstrap approach has a similar shape to that using the normal distribution but a higher mean value. Based on these observations from the top panels, we expect the price of a future written on $\text{CWSI}_{[t_0,t_1]}^{(1)}$ not to differ significantly when the distribution is switched from normal to HYP. However, the price of an option with a certain strike price may experience a larger price difference since the tails of the density curves under the two distributions are significantly different. The semi-parametric model is not likely to reduce pricing errors arising from misspecified distribution for futures written on $\text{CWSI}_{[t_0,t_1]}^{(1)}$ .

The bottom panels demonstrate that the mean values of simulated $\text{CWSI}_{[t_0,t_1]}^{(2)}$ under the normal assumption are much higher than those under the HYP assumption. In Fig. 4, we observe that the normal distribution leads to higher densities for wind speed between 3 and 4.88 and slightly lower densities for wind speed above 4.88 in comparison to the HYP distribution. Since $\text{CWSI}_{[t_0,t_1]}^{(2)}$ only cumulates wind speed between 3 and 25 and the probability of wind speed above 4.88 is very low, the normal assumption results in significantly higher mean values for $\text{CWSI}_{[t_0,t_1]}^{(2)}$ . The density curve using the bootstrap approach has a shape similar to that using the normal distribution but a slightly lower mean value. Based on the observations made on the bottom panels, we expect large price differences for derivatives written on $\text{CWSI}_{[t_0,t_1]}^{(2)}$ under the normal and hyperbolic assumptions. The price difference between using the bootstrap and HYP assumption is expected to be smaller. Therefore, the semi-parametric model may reduce pricing errors arising from misspecified distribution by a small amount.

7.2 Pricing results

The tick size $D$ is set to 1 without loss of generality. The 25th, 50th, and 75th percentiles of the historical wind speed index values in 2012–2021 are used as the strike prices for put options. The risk-free interest rate $r$ is set at 4% per annum. Tables 5 and 6 present the prices of derivatives written on $\text{CWSI}_{[t_0,t_1]}^{(1)}$ and $\text{CWSI}_{[t_0,t_1]}^{(2)}$ respectively. In both tables, we show the prices obtained with the two proposed pricing approaches and under various choices of distributions and $\theta$ values. When $\theta =0$ , the two probability measures $P$ and $Q$ are the same, and thus the two pricing approaches yield exactly the same prices. When $\theta \lt 0$ , the two pricing approaches result in fairly small differences under the same distribution assumption, thereby validating the use of both approaches. As $\theta$ becomes more negative, the higher risk premiums demanded by investors lead to lower future prices and higher put option prices.

Table 5. Prices of futures and options written on $\text{CWSI}_{[t_0,t_1]}^{(1)}$ under various choices of distributions and $\theta$ values

Table 6. Prices of futures and options written on $\text{CWSI}_{[t_0,t_1]}^{(2)}$ under various choices of distributions and $\theta$ values

In Table 5, the differences of prices obtained under the normal and HYP assumptions are very small for the future contract and put options with $K=236.11$ and $K=229.62$ . The differences in prices obtained using the bootstrap approach and the HYP distribution are larger for these derivatives. This finding is consistent with the observation that the mean values of simulated $\text{CWSI}_{[t_0,t_1]}^{(1)}$ are very close under the two parametric distributions but higher using bootstrap from Fig. 5.

For the put options with $K=201.40$ , the HYP assumption and the bootstrap yield similar prices, with both being lower than the prices obtained under the normal assumption. In Fig. 5, we observe that the densities of simulated $\text{CWSI}_{[t_0,t_1]}^{(1)}$ for $\text{CWSI}_{[t_0,t_1]}^{(1)}\lt 201.40$ are generally higher under the normal assumption than they are under the HYP assumption. Therefore, the probability of non-zero option payoffs is also higher under the normal assumption, resulting in a higher option price. On the other hand, the left tails of the density curves under the HYP assumption and the bootstrap approach are close to each other, thereby explaining the close prices under the two assumptions. Therefore, the left tails of the density curves in Fig. 5 play a key role in the price of this option.

Table 5 indicates that the pricing error arising from misspecified distribution (normal) is small for futures and put options with relatively high strikes written on $\text{CWSI}_{[t_0,t_1]}^{(1)}$ . The semi-parametric model can only reduce the pricing error for options with low strikes.

In Table 6, the normal and HYP assumptions lead to very large price differences for all the derivatives written on $\text{CWSI}_{[t_0,t_1]}^{(2)}$ . Fig. 5 shows that the mean values of simulated $\text{CWSI}_{[t_0,t_1]}^{(2)}$ under the normal assumption are much higher than those under the HYP distribution. As a result, the normal assumption leads to higher prices for futures and lower prices for put options. The differences in prices obtained under the HYP assumption and the bootstrap approach are smaller than those obtained under the normal and HYP distributions. Therefore, the pricing error arising from misspecified distribution is very high for derivatives written on $\text{CWSI}_{[t_0,t_1]}^{(2)}$ . The semi-parametric model can partially reduce the pricing error.

8. Conclusion

In this paper, we investigate how model risk, in particular misspecification of model distribution, affects wind derivative prices. We first estimate s-AR-s-GARCH models with four different parametric distributions – normal, NIG, HYP, and GHYP. The NIG, HYP, and GHYP distributions are found to yield a much better fit and significantly different parameter estimates than the normal distribution. Therefore, the bias of the QML estimator, which maximizes normal likelihood, is high when fitting wind speed data.

We develop two simulation approaches based on the conditional Esscher transform to determine the risk-neutral prices of wind speed derivatives. Assuming normal, hyperbolic, and non-parametric distributions for wind speed, the prices of wind speed futures and put options written on two different wind speed indices are calculated. The price difference between using the normal and hyperbolic distributions can be used as a measure of model risk. Our pricing results show that the extent of model risk depends on the choice of wind speed index and the derivative payoff structure.

Significant price differences are observed when the derivative is written on the wind speed index that approximates wind power production, regardless of the derivative payoff structure. The semi-parametric wind speed model, which uses the bootstrap procedure to incorporate non-normality, can only help reduce pricing error by a small amount. Since the semi-parametric wind speed model also uses the QML estimates, we can conclude that the pricing errors predominantly arise from the bias in the QML estimates. Minor pricing differences have resulted when the wind speed index is a simple sum of the daily average wind speed. However, the differences are more notable for put options with low strike prices.

Without considering other weather conditions that may affect wind power production, the NAREX WIDE index can be viewed as a function of wind speed. The pricing approaches developed in this paper can be easily adapted to price wind power futures written on the NAREX WIDE. Since the NAREX WIDE index is determined by wind conditions and power curve, it is similar to the wind speed index $\text{CWSI}_{[t_0,t_1]}^{(2)}$ considered in this paper. Since we have observed high model risk in the pricing of derivatives written on $\text{CWSI}_{[t_0,t_1]}^{(2)}$ , the model risk is expected to be high in the pricing of wind power futures as well. Therefore, choosing a suitable distribution is essential for accurately pricing wind power futures.

It is important to note that the wind derivatives constructed in this paper are hypothetical and intended for illustrative purposes; thus, no actual traded prices can be observed. In future research, we plan to apply the proposed pricing approaches to the traded wind power futures written on the NAREX WIDE index and evaluate the performance of the proposed pricing method and the risk of model misspecification.

Actuarial expertise holds significant value in the quantification and management of risks inherent in renewable energy markets, such as wind power. In this paper, we demonstrate how actuaries can employ their statistical modeling skills to refine wind speed models, as well as leverage their financial analysis expertise to devise pricing methodologies for related financial instruments. Furthermore, actuaries possess the necessary skills to assess and manage risks associated with these financial instruments. There are several potential avenues for future research in the renewable energy market that can greatly benefit from actuarial expertise.

Firstly, future research should focus on assessing the hedge effectiveness of wind power futures. The NAREX WIDE index reflects the overall power output in Germany, while smaller energy producers operate wind turbines in specific geographical locations. Basis risk arises from the mismatch between national wind power output and output from individual energy producers, potentially reducing the effectiveness of wind power futures. It is crucial to quantify the extent to which existing wind power futures can protect against variability in wind power production.

Secondly, future research should aim to develop optimal hedging strategies for energy producers utilizing wind derivatives. This will enable them to better manage revenue uncertainties and safeguard against the financial impacts of unfavorable wind conditions.

Lastly, actuaries should collaborate with professionals from diverse fields, including climate scientists, economists, and energy experts, to foster a comprehensive understanding of the challenges and opportunities linked to renewable energy and its associated financial risks. This interdisciplinary approach can lead to the creation of innovative financial instruments and risk management strategies tailored to the specific needs of the renewable energy sector.

Integrating actuarial expertise into risk quantification and management within the renewable energy sector will lead to a more resilient and efficient market for clean energy. This, in turn, will bolster the global transition towards sustainable energy production.