Journal of Energy Markets

Risk.net

Optimal management of green certificates in the Swedish–Norwegian market

Fred Espen Benth, Marcus Eriksson and Sjur Westgaard

  • Stochastic model of green certificate prices in the Swedish-Norwegian market.
  • Optimal management of green certificates for a producer of renewable power.
  • Stochastic optimal control and dynamic programming with singular controls.

We propose and investigate a valuation model for the income of selling tradeable green certificates (TGCs) in the Swedish–Norwegian market, formulated as a singular stochastic control problem. Our model takes into account the production rate of renewable energy from a “typical” plant, the price of TGCs and the cumulative amount of certificates sold. We assume that the production rate has a dynamics given by an exponential Ornstein–Uhlenbeck process, and the logarithmic TGC price has a dynamics given by a Lévy process. For this class of dynamics, we find optimal decision rules for the state variables and a closed-form solution to the control problem. A case study of ICAP prices and wind production data from Denmark backs up our model choice and shows the relevance of this pricing approach.

1 Introduction

The need for renewable energy has increased a lot in the last decade as the awareness of climate change has increased. The renewable energy technologies are often very expensive and require external financial support to be realized. An alternative to direct government funding is the development of a market for green certificates, also called tradeable green certificates (TGCs). The certificates are purely financial objects used to reach a desirable production capacity of electricity from renewable resources. The idea is that the end consumers finance the renewable energy technologies by purchasing TGCs.

The producer of “green” electricity has the right to sell one certificate per unit produced, while the retail electricity providers (eg, the consumers in this market) are obliged to cover their share of electricity consumption from renewable sources. This is done by purchasing green certificates. The producer is given a number of certificates by the government, based on expected production, to sell in the market. In this way a market-based subsidy is introduced, creating a direct link between power consumption and “green” energy generation.

The Swedish electricity certificate market was established in 2003, and Norway enrolled in the market in 2012. Before the common market between Sweden and Norway, 13.3 TWh was financed via the Swedish certificate market. The common goal was to increase renewable electricity production by 26.4 TWh between 2012 and 2020, shared equally between Sweden and Norway. The market will continue until 2035, with a target of completely financing 198 TWh of new renewable electricity production (NVE 2013). For comparison, Norway and Sweden produced a total of 282 TWh in 2014.

The aim of this paper is to find the value of an income stream obtained by selling certificates in the particular case of the Swedish–Norwegian green certificate market. By observing the market price of the certificates, the holder can maximize their income stream by finding the best amount to sell at the best time. We formulate the optimization problem as a finite-horizon singular stochastic control problem, with the state variables being the certificate spot price, the production rate of power from renewable sources and the number of certificates sold. We allow certificates to be sold in small quantities over time, or all at once, which includes singular stochastic controls in our optimization problem.

Since one of the aims of TGCs is to cover the producer’s losses due to the high installation costs of renewable energy technologies, the consumer pays a higher price for electricity, namely the sum of the market prices for electricity and TGCs. Both the electricity price and the TGC price will be driven by the demand for power, and we can expect a dependency between the two.

In stochastic modeling of financial prices, exponential Lévy processes are commonly used (see Cont and Tankov 2004; Schoutens 2003). Such processes provide a great flexibility in capturing the stylized features of financial price data, such as heavy tails, skewness and kurtosis in the distribution of prices. To obtain a flexible model for the TGC price, we will assume that it follows an exponential Lévy process, where logarithmic returns are modeled by the normal inverse Gaussian (NIG) distribution (see Barndorff-Nielsen (1998) for a financial application of the NIG distribution). As it turns out, this model will fit observed TGC price data rather well, as we will see in an empirical case study. We mention in passing that Lévy processes are frequently used in modeling electricity prices (see Benth et al 2008).

The production rate is highly influenced by weather factors such as rain, sunshine and wind, given renewable power generators such as hydroelectric power, photovoltaic cells and wind turbines. It is reasonable to assume that these weather factors are stationary, varying around some seasonal mean level. To describe the statistical features of wind speed dynamics in discrete time, the most commonly used models are autoregressive moving-average (ARMA) time series models. These have also been used for modeling time series for temperature. The analog for continuous time is the continuous-time ARMA (CARMA) process. Applications of these, and other related weather models, in the weather market can be found in Benth and Saltyte Benth (2013) and the references therein. In particular, Benth and Saltyte Benth suggest a CARMA model for temperature and an exponential CARMA model for wind speed. Motivated by this, we will assume that the production rate follows an exponential Ornstein–Uhlenbeck (OU) process. The OU process is stationary, and a special case of the more general CARMA processes. However, it is also Markovian, which is theoretically convenient when analyzing our optimization problem.

Previous work relating to green certificates markets is rather scarce. The term paper by Goldstein (2010) gives an overview of the market and its implications from a political and economical point of view, as well as a discussion on the equilibrium price and the Swedish–Norwegian market. There is a stream of literature on price determination in the TGC market, including fundamental equilibrium and agent-based models (see, for example, Amundsen and Mortensen 2001; Aune et al 2012; Unger and Ahlgren 2005; Wolfgang et al 2015). The impact of regulatory changes on the price volatility of TGC is analyzed in Fagiani and Hakvoort (2014). Frogner and Hustveit (2015) study green investment decisions using dynamic programming. The related solar renewable energy certificates (SREC) market in the United States is investigated by Coulon et al (2015), who focus on understanding the price dynamics and propose a structural model for renewable energy certificates that is able to incorporate important features.

To the best of our knowledge, there are no papers discussing how a producer should manage TGCs optimally in a continuous-time market context. We provide a general framework for a valuation model for selling TGCs optimally, as well as a general model for the underlying TGC price dynamics. A closed-form solution to the valuation model is provided and explicitly calculated when the logarithmic TGC price process is NIG distributed. We also conduct an empirical analysis, which demonstrates that our proposed model and the NIG distribution are very suitable for work on price dynamics. In this case, we also calculate the numerical value of the contract based on the empirical data.

We also highlight the flexibility of the valuation model, as it can be used for any exponential Lévy process with finite moments to model the TGC price dynamic, yet remains analytically tractable.

The paper is organized as follows. In Section 2, we give the framework for our valuation model, and introduce the dynamics of the production rate. In Section 3, we introduce the exponential Lévy process for the TGC price dynamics and derive a Hamilton–Jacobi–Bellman (HJB) equation for the valuation problem. We then derive criteria for the optimal strategy, and show the optimality via a verification theorem. The main result of the paper gives a rather explicit solution to the valuation model. In Section 4, we conduct a case study followed by an explicit numerical calculation of the optimal value. For the numerical part we conduct an empirical analysis of TGC spot price data. Proofs and intermediate results are collected in Section 5. We offer some conclusions in Section 6.

2 The singular stochastic control problem

Let (Ω,,𝔽,𝑷), where 𝔽={t}t0, be a complete filtered probability space satisfying the usual conditions. We denote by t:={𝑿(u),ut} the σ-algebra generated by the state process 𝑿(u). Also, we assume that the state process has the Markov property.

We formulate as follows the singular stochastic control problem for optimal management of the green certificates held by the producer. Let X(t) be the price of a green certificate at time t, and denote by P(t) the accumulated production of “clean power” from a producer entitled to receive certificates. We introduce R(t) as the production rate at time t, and thus

  dP(t)=R(t)dt.  

Furthermore, we denote by A(t) the cumulative number of certificates sold up to time t.

The market is organized such that the producer of “clean power” is granted a number of certificates proportional to the production on a regular basis. We approximate this as a continuous income of certificates proportional to the production rate R(t). Hence, at time t we have cP(t) accumulated certificates obtained from production, with c>0 being the proportionality constant. These certificates can, once received, be sold at any later time. Our selling strategy is modeled by the control A(t). We assume the following conditions hold on the set of controls.

  • (C)

    A is a positive, nondecreasing and adapted stochastic process, with paths being right continuous with left limits. We let A(0-)=0. In addition, A(t)cP(t) for all tT. We call such controls “admissible”, and denote the set of admissible As by 𝒜.

Note that T is a finite trading horizon, typically the total duration of the certificate market. The condition A(t)cP(t) prohibits short selling of certificates.

We introduce the process Z(t), which measures the number of certificates held at time t, ie,

  Z(t)=cP(t)-A(t).  

We observe that, for any A𝒜, it holds that Z(t)0 for all tT.

With the above notation, we have the state variable

  𝑿(t)=(X(t),R(t),Z(t))  

controlled by A𝒜. The expected value of the income flow from selling certificates becomes

  J(t,x,ϱ,z:A)=𝔼[tTe-r(s-t)X(s)dA(s)|𝑿(t-)=(x,ϱ,z)]  

for any A𝒜(t), where 𝒜(t) is the set of admissible controls and time starts at t. We have denoted by r>0 the constant discount rate. Note that, as A is monotonically nondecreasing, it is of finite variation on the interval [t,T]. Hence, the integral with respect to A inside the expectation operator above is interpreted in the Lebesgue–Stieltjes sense. Our stochastic control problem is now to find an optimal A^𝒜(t) such that

  V(t,x,ϱ,z):=supA𝒜(t)J(t,x,ϱ,z:A)=J(t,x,ϱ,z:A^).   (2.1)

We analyze this singular stochastic control problem by the method of dynamic programming.

We observe that if t=T, the optimal control is to sell all the certificates that the producer holds. Hence, if Z(T-)=z, the optimal control is ΔA^(T)=z. The value of selling these certificates is then given by

  V(T,x,ϱ,z)=𝔼[X(T)ΔA^(T)𝑿(T-)=(x,ϱ,z)]=xz.   (2.2)

This provides us with a terminal condition for the value function.

Note that we set the optimization problem under the market probability P rather than taking the expectation with respect to any risk-neutral probability Q. Indeed, in mathematical finance, the reason for pricing derivatives by taking the expected payoff with respect to a risk-neutral probability QP comes from the fact that we can hedge the derivatives by the underlying asset. In our context, the producers have been granted American-type options, where the underlyings are the TGC and production. The production is an external variable that cannot be used for hedging in a financial sense (we cannot “trade the wind”), whereas the spot market for TGCs is far from functioning as a liquid financial market. The retailers purchase TGCs by need, and not for speculative purposes. Furthermore, in our setup there is a short-selling constraint that creates additional incompleteness in the market. Taking all these aspects into consideration, we are in reality in a highly incomplete market. We have therefore chosen to state the optimization problem under the market probability P, where r is a discount factor not necessarily equal to the so-called risk-free interest rate. However, we emphasize that it is rather simple to modify our analysis to accommodate a liquid market situation by assuming a drift in the TGC price model (see, for example, Cont and Tankov 2004). The change in the dynamics for the production rate under such a change of probability is, on the other hand, still questionable even in a liquid TGC market.

We will focus our optimal control problem on some particular classes of state processes X and R of practical relevance and interest. As the production rate R is highly influenced by weather factors, which are stationary, varying around some seasonal mean level, a simple, yet natural model is to assume that the dynamics of R follows an exponential OU process:

  R(s)=eU(s),   (2.3)

where U(s) is a mean-reverting OU process driven by a Brownian motion:

  dU(s)=(μ-αU(s))ds+σudBu(s),U(t)=ln(R(t)).  

Then the dynamics of R(s) reads as

  dR(s)=aR(R(s))R(s)dt+σuR(s)dBu(s),R(t)=ϱ,  

where

  aR(R(s)):=μ-αln(R(s))+12σu2.   (2.4)

The constants μ, α and σu denote the mean-reversion level, the rate of mean reversion and the volatility of the process U(s). Bu is a Brownian motion, where the superscript u indicates that it is related to the process U(s). The explicit solution to R(s), starting at time t, is given by

  R(s)=exp(ln(R(t))e-α(s-t)+μα(1-e-α(s-t))+tsσue-α(s-v)dBu(v)).   (2.5)

In the next section, we specify the price process, X.

3 A price model for TGC and dynamic programming

In this section, we find a closed-form solution for the optimal value (2.1). We assume the logarithmic price, denoted by Y, is a Lévy process with finite moments. The production rate R(s) is assumed to be an exponential OU process given by (2.5). As a by-product, we also obtain the optimal strategy, ie, the optimal control. The result is concluded in Theorem 3.8. We start by specifying the price model for X.

3.1 The TGC price model

Let

  X(s)=xexp(Y(s)),  

where the dynamics of the Lévy process Y is given by

  dY(v)=γdv+σYdBY(v)+|ξ|<1ξN~(dv,dξ)+|ξ|1ξN(dv,dξ),   (3.1)

and BY is a Brownian motion correlated with Bu with correlation coefficient ρ. N is a Poisson random measure with Lévy measure ν(dξ) as compensator. Furthermore, assume that X has finite moments, ie, we suppose that the condition

  |ξ|1ek|ξ|ν(dξ)<   (3.2)

holds for some k>2. As a consequence, we have

  {0}|eξ-1-ξ|ν(dξ)<   (3.3)

and

  |ξ|1|ξ|ν(dξ)<.   (3.4)

We can write

  |ξ|<1ξN~(dv,dξ)+|ξ|1ξN(dv,dξ)=|ξ|<1ξN~(dv,dξ)+|ξ|1ξN~(dv,dξ)+|ξ|1ξν(dξ).  

Hence,

  dY(v)=γ~dv+σYdBY(v)+{0}ξN~(dv,dξ),  

where

  γ~:=γ+|ξ|1ξν(dξ).  

By the Ito formula for semimartingales, we obtain the dynamics for X(v):

  dX(v)=aXX(v)dv+σYX(v)dBY(v)+{0}X(v-)(eξ-1)N~(dv,dξ),  

where

  aX:=γ~+12σY2+{0}(eξ-1-ξ)ν(dξ).   (3.5)
Lemma 3.1.

Let Y(1) have Lévy triplet (γ,σY2,ν(dξ)) and characteristic function ϕ(u). Suppose that (3.3) holds. Then,

  aX=lnϕ(-i),  

where i is the imaginary unit.

Proof.

We have that

  aX:=γ+|ξ|1ξν(dξ)+{0}(eξ-1-ξ)ν(dξ)+12σY2.  

By the Lévy–Khintchine formula for the logarithm of the characteristic function, we have

  lnϕ(u) =iγu-12σY2u2+{0}(eiuξ-1-iuξ𝟏(|ξ|<1))ν(dξ)  
    =iγu-12σY2u2+{0}(eiuξ-1-iuξ𝟏(|ξ|<1)+iuξ𝟏(|ξ|>1)-iuξ𝟏(|ξ|>1))ν(dξ)  
    =iγ~u-12σY2u2+{0}(eiuξ-1-iuξ)ν(dξ)  

for u. Here, we used (3.4) in the second line. By condition (3.3), u can be extended to the complex plane. Taking u=-i yields

  lnϕ(-i)=γ~+12σY2+{0}(eξ-1-ξ)ν(dξ).  

Hence, the result follows. ∎

3.2 The valuation model

We will now solve the control problem defined in (2.1). First, we define the space (t,T,BY,Bu) as all F(t,x,ϱ,z)C1,2,2,1 such that, for any admissible control A𝒜(t), the processes

  θ tθe-rsFx(s,X(s),R(s),Z(s))σ1X(s)dBY(s),  
  θ tθe-rsFϱ(s,X(s),R(s),Z(s))σ2R(s)dBu(s)  

are martingales, for tθT. The HJB equation is now derived via Bellman’s principle of optimality.

Proposition 3.2.

Suppose that V(t,x,ϱ,z)(t,T,BY,Bu), and that the process

  θtθ{0}e-rs[V(s,X(s)eξ,R(s),Z(s))-V(s,X(s),R(s),Z(s))]N~(ds,dξ)  

is a martingale. Then, for all t[0,T], the corresponding HJB equation associated with the value function V is

  max(Vt+V-rV,-Vz+x)=0,   (3.6)

where the operator , acting on functions FF(t,x,ϱ,z)C1,2,2,1, is defined as

  F :=aXxFx+aR(ϱ)ϱFϱ+cϱFz+12σY2x2Fxx+12σu2ϱ2Fϱϱ+ρσYσuxϱFxϱ  
      +{0}[F(t,xeξ,ϱ,z)-F(t,x,ϱ,z)-x(eξ-1)Fx(t,x,ϱ,z)]ν(dξ),   (3.7)

with aR(ϱ) given by (2.4), and aX by (3.5).

Proof.

See Section 5. ∎

To proceed, the following lemma will be useful for further calculations.

Lemma 3.3.

Let ϕ(u) be the characteristic function of Y(1), defined in (3.1), and suppose the condition of Lemma 3.1 holds. Then, for all s[t,T], the following expected equalities hold.

  1. (i)

    We have

      𝔼[X(s)R(s)𝑿(t)=(x,ϱ,z)]=xexp(ln(ϱ)e-α(s-t)+μα(1-e-α(s-t))+σu24α(1-e-2α(s-t))(2-ρ2)+σuσYα(1-e-α(s-t))+aX(s-t)).   (3.8)
  2. (ii)

    We have

      𝔼[X(T)𝑿(t)=(x,ϱ,z)]=xexp[aX(T-t)].   (3.9)
  3. (iii)

    We have

      𝔼[X(T)R(s)𝑿(t)=(x,ϱ,z)]=xexp(ln(ϱ)e-α(s-t)+μα(1-e-α(s-t))+σu24α(1-e-2α(s-t))(2-ρu2)+σuσYα(1-e-α(s-t))+aX(T-t)).   (3.10)
Proof.

See Section 5. ∎

Define

  h h(t,ϱ,s):=1x𝔼[X(s)R(s)𝑿(t)=(x,ϱ,z)],  
  h~ h~(t,ϱ,s):=1x𝔼[X(T)R(s)𝑿(t)=(x,ϱ,z)],  

and

  H H(t,ϱ,T):=tTe-r(s-t)h(t,ϱ,s)ds,   (3.11)
  H~ H~(t,ϱ,T):=tTe-r(s-t)h~(t,ϱ,s)ds.   (3.12)

Consider the admissible control A~1, defined by ΔA~1(t)=z and dA~1(s)=cR(s)ds for s>t. Define

  Φ(t,x,ϱ,z) :=J(t,x,ϱ,z:A~1)  
    =𝔼[tTe-r(s-t)X(s)dA~1(s)|𝑿(t-)=(x,ϱ,z)].   (3.13)

Then,

  Φ(t,x,ϱ,z) =𝔼[xz+tTe-r(s-t)X(s)cR(s)ds|𝑿(t)=(x,ϱ,z)]  
    =xz+ctTe-r(s-t)𝔼[X(s)R(s)𝑿(t)=(x,ϱ,z)]ds  
    =xz+cxH(t,ϱ,T).  
Proposition 3.4.

If

  aXr,   (3.14)

then Φ(t,x,ϱ,z) defined in (3.2) solves the HJB equation (3.6). Furthermore, Φ(t,x,ϱ,z) is dominated by the value function, ie,

  Φ(t,x,ϱ,z)V(t,x,ϱ,z),  

and Φ(T,x,ϱ,z)=V(T,x,ϱ,z)=xz.

Proof.

From the relation

  Φ(t,x,ϱ,z)=xz+cxH(t,ϱ,T),  

it is clear that Φ(T,x,ϱ,z)=xz. Also, we see that x-Φz=0. It remains to show that Φt+Φ-rΦ0. Using the definition of H(t,ϱ,T) in (3.11), we calculate

  Φt+Φ-rΦ =cxtTe-r(s-t)(rh+ht)ds-cxϱ+aXx(z+ctTe-r(s-t)hds)  
      +cxaR(ϱ)ϱtTe-r(s-t)hϱds+cxσYσuρϱtTe-r(s-t)hϱds  
      +cxϱ+12cxσu2ϱ2tTe-r(s-t)hϱϱds-r(xz+cxtTe-r(s-t)hds)  
      +{0}[xeξz+cxeξH(t,ϱ,T)-(xz+cxH(t,ϱ,T))-x(eξ-1)(z+cH(t,ϱ,T))]ν(dξ)  
    =(aX-r)xz+cxtTe-r(s-t)[ht+aXh+aR(ϱ)ϱhϱ+ρσYσuϱhϱ+12σu2ϱ2hϱϱ]ds.  

Inserting the derivatives of h, given in Section 5, yields

  Φt+Φ-rΦ =(aX-r)xz+cxtTe-r(s-t)h[M(s,t)+aX+e-α(s-t)(aR(ϱ)+ρσYσu-12σu2)+12σu2e-2α(s-t)]ds.  

The first term is clearly nonpositive due to (3.14), since xz>0. For the integrand, from the expressions for aR(ϱ)=μ-αln(ϱ)+12σu2 and M(s,t) in (5.9), we obtain

    M(s,t)+aX+e-α(s-t)(aR(ϱ)+ρσYσu-12σu2)+12σu2e-2α(s-t)  
            =e-α(s-t)[αln(ϱ)-μ-σuσY]+12σu2(ρ2-2)e-2α(s-t)-12σY2-lnϕ(-i)  
              +aX+e-α(s-t)(μ-αln(ϱ)+12σu2+ρσYσu-12σu2)+12σu2e-2α(s-t)  
            =e-α(s-t)[σuσY(ρ-1)+12σu2(ρ2-2)e-α(s-t)]+aX-aX  
            0.  

The inequality follows since σuσY(ρ-1)+12σu2(ρ2-2)e-α(s-t)0. Hence, Φt+Φ-rΦ0 since cx0. The domination follows since A~1 is admissible. ∎

If aX>r, the control A~1 would violate the HJB equation. Consider instead the admissible control A~2 defined as A~2(s)=0 for s[t,T) and

  ΔA~2(T)=Z(T-)=z+tTcR(s)ds.  

Set

  Φ(t,x,ϱ,z)=J(t,x,ϱ,z:A~2).   (3.15)

Then,

  Φ(t,x,ϱ,z) =𝔼[e-r(T-t)X(T)ΔA~(T)𝑿(t-)=(x,ϱ,z)]  
    =e-r(T-t)𝔼[X(T)(z+tTcR(s)ds)|𝑿(t)=(x,ϱ,z)]  
    =e-r(T-t)z𝔼[X(T)𝑿(t)=(x,ϱ,z)]  
      +e-r(T-t)ctT𝔼[X(T)R(s)𝑿(t)=(x,ϱ,z)]ds.  

Thus, by Lemma 3.3(ii) and (3.12), we obtain

  Φ(t,x,ϱ,z)=xze(aX-r)(T-t)+cxH~(t,ϱ,T).   (3.16)
Proposition 3.5.

If

  raX,   (3.17)

then Φ(t,x,ϱ,z) defined in (3.15) solves the HJB equation (3.6). Furthermore, Φ(t,x,ϱ,z) is dominated by the value function, ie,

  Φ(t,x,ϱ,z)V(t,x,ϱ,z),  

and Φ(T,x,ϱ,z)=V(T,x,ϱ,z)=xz.

Proof.

By (3.16), it is clear that the terminal condition holds. Turning to the HJB equation, we have

  x-Φz=x(1-e(aX-r)(T-t))0   (3.18)

by condition (3.17). Similarly to the proof of Proposition 3.4, we get

  Φt+Φ-rΦ =[e(aX-r)(T-t)(r-aX)]xz+cxtTe-r(T-t)[h~t+aXh~+aR(ϱ)ϱh~ϱ+ρσYσuϱh~ϱ+12σu2ϱ2h~ϱϱ]ds.  

The first term is nonpositive by (3.17). The second is also nonpositive by similar calculations to those in the proof of Proposition 3.4. The derivatives of Φ and h~ can be found in Section 5. The domination follows, since A~2 is an admissible control. ∎

Theorem 3.6.

Suppose that Φ(t,x,ϱ,z)(t,T,BY,Bu), and that the process

  θtθ{0}e-rs[Φ(s,X(s)eξ,R(s),Z(s))-Φ(s,X(s),R(s),Z(s))]N~(ds,dξ)   (3.19)

is a martingale. If Φ solves the HJB equation (3.6), with Φ(T,X(T),R(T),Z(T))=V(T,X(T),R(T),Z(T)), then ΦV for all (t,x,ϱ,z)[0,T]×+×+×[0,M].

Proof.

For tθT, by Ito’s formula we obtain

  tθd(e-rsΦ(s,X(s),R(s),Z(s))) =tθe-rsΦt(s,X(s),R(s),Z(s))-re-rsΦt(s,X(s),R(s),Z(s))ds  
      +tθe-rsΦx(s,X(s),R(s),Z(s))aXX(s)ds  
      +tθe-rsΦϱ(s,X(s),R(s),Z(s))aR(R(s))R(s)ds  
      +tθe-rsΦz(s,X(s),R(s),Z(s))cR(s)ds  
      +tθe-rsΦx(s,X(s),R(s),Z(s))σYX(s)dBY(s)  
      +tθe-rsΦϱ(s,X(s),R(s),Z(s))σuR(s)dBu(s)  
      -tθe-rsΦz(s,X(s),R(s),Z(s))dA(s)  
      +tθe-rsΦxx(s,X(s),R(s),Z(s))12σY2X2(s)ds  
      +tθe-rsΦϱϱ(s,X(s),R(s),Z(s))12σu2R2(s)ds  
      +tθe-rsΦxϱ(s,X(s),R(s),Z(s))σYσuρX(s)R(s)ds  
      +tθ{0}e-rs[Φ(s,X(s-)eξ,R(s-),Z(s-))-Φ(s,X(s-),R(s-),Z(s-))]N~(ds,dξ)  
      +tθ{0}e-rs[Φ(s,X(s)+X(s-)(eξ-1),R(s),Z(s))-Φ(s,X(s),R(s),Z(s))-X(s)(eξ-1)]ν(dξ)ds.  

Taking the expectation and using the conditions in the theorem yields

    Φ(t,x,ϱ,z)-e-r(θ-t)𝔼[Φ(θ,X(θ),R(θ),Z(θ))𝑿(t)=(x,ϱ,z)]  
        =𝔼[tθe-r(s-t)Φz(s,X(s),R(s),Z(s))dA(s)+tθe-r(s-t)[-(Φt+Φ-rΦ)(s,X(s),R(s),Z(s))]ds|𝑿(t)=(x,ϱ,z)].  

Take θ=T. Then,

  Φ(t,x,ϱ,z) =𝔼[tTe-r(s-t)Φz(s,X(s),R(s),Z(s))dA(s)  
          +tTe-r(s-t)[-(Φt+Φ-rΦ)(s,X(s),R(s),Z(s))]ds  
              +e-r(T-t)Φ(T,X(T),R(T),Z(T))|𝑿(t)=(x,ϱ,z)].  

Since Φ satisfies the HJB equation, at time s we have

  Φz(s,X(s),R(s),Z(s))X(s),-(Φt+Φ-rΦ)(s,X(s),R(s),Z(s))0,  

and by assumption we have

  Φ(T,X(T),R(T),Z(T))=V(T,X(T),R(T),Z(T))0.  

It follows that

  Φ(t,x,ϱ,z)𝔼[tTe-r(s-t)X(s)dA(s)|𝑿(t)=(x,ϱ,z)].  

Since this inequality holds for any admissible control, it also holds for the supremum over all such controls. Hence,

  Φ(t,x,ϱ,z)V(t,x,ϱ,z).  

Lemma 3.7.

For the functions Φ(t,x,ϱ,z) defined in (3.2) and (3.15), condition (3.19) holds and Φ(t,x,ϱ,z)(t,T,BY,Bu).

Proof.

See Section 5. ∎

From the Verification Theorem (Theorem 3.6), we conclude the following.

Theorem 3.8.

For aX<r,

  V(t,x,ϱ,z)=xz+cxH(t,ϱ,T),   (3.20)

and the optimal control is A~1.

For raX,

  V(t,x,ϱ,z)=xze(aX-r)(T-t)+cxH~(t,ϱ,T),   (3.21)

and the optimal control is A~2.

Proof.

Since A~1 and A~2 are admissible, the conclusion follows directly from Propositions 3.4 and 3.5 together with Lemma 3.7 and Theorem 3.6. ∎

The optimal value of the income from selling certificates is thus given by Theorem 3.8, and as a by-product we get the optimal strategy. Theorem 3.8 also provides us with conditions on which strategy to use, depending on the sign of aX-r, ie, A~1 or A~2. Hence, in practice we only need to determine the sign of aX-r. Then, by Theorem 3.8, a solution for the optimal value is given by (3.20) or (3.21). Apart from an integration in definitions of H and H~, we have a closed-form solution.

Furthermore, note the interpretation of the optimal strategies A~1 and A~2: if the expected rate of return is less than the discount rate r, the TGC holder should sell all certificates they hold at time t and continue selling at the same rate as certificates are obtained from the production of renewable energy. On the other hand, if the expected rate of return is greater than r, the TGC holder should wait, and then sell all certificates at maturity.

In the next section, we will conduct a case study where we explicitly calculate the optimal value. First, just to illustrate the simplicity of what control to choose, consider the following simple example.

Figure 1: Time series for traded TGC spot price data X
Figure 1: Time series for traded TGC spot price data X. Data collected from February 17, 2009 to November 1, 2016.
Example 3.9.

Let the logarithmic price follow a Brownian diffusion process:

  dY(t)=γdt+σYdBY(t).  

This results in a price X having geometric Brownian motion dynamics

  dX(v)=aXX(v)dv+σYX(v)dBY(v),  

where aX=γ+12σY2. Consider a time series of certificate spot prices collected from ICAP.11The data is provided by Montel. In Figure 1, we show the time series of traded prices from February 17, 2009 to November 1, 2016. All together, there are 401 prices denoted in SEK/MWh. The certificates are not traded every day, and in the time span considered here we have trades on average about every fourth day when the exchange is open, or five trades per month. There is a general downward trend in the prices, where the annual value of the drift aX is estimated to be -8.5% and the sign of aX-r is therefore naturally negative. Then, the optimal strategy is found to be A~1 by Theorem 3.8.

4 A case study

In this section, we perform a case study that illustrates our theoretical analysis of the optimal management of certificates using actual prices and production data. We base our investigations on the ICAP price data introduced above, and on realized wind power production from Denmark. The latter is chosen in order to obtain reasonable estimates on the model for the production rate, as we do not have available production data from specific wind power plants in Norway or Sweden.

Figure 2: Empirical density of the ICAP log returns
Figure 2: Empirical density of the ICAP log returns (solid line) with the MLE fitted normal (dotted) and NIG (dashed line) distributions.
Table 1: MLE estimated NIG parameters.
Parameter MLE
α^NIG 22.87
β^NIG 0-5.18
μ^NIG 0.00065
δ^NIG 0.0047

Let us first have a closer look at the TGC prices. In the example at the end of the previous section, we discussed a geometric Brownian motion model. However, in Figure 2 we show the empirical density of the log returns together with the fitted normal distribution on a logarithmic frequency scale. As is evident, the tails are not explained well at all, but, using the more flexible NIG distribution (see, for example, Barndorff-Nielsen (1998) for the NIG distribution applied in finance), we obtain a much better fit.

The NIG distribution is a four-parameter family of distributions, with parameters denoted αNIG, βNIG, δNIG and μNIG. The location of the distribution is expressed by μNIG, and the scale by δNIG. The skewness parameter is βNIG, and the tail heaviness is modeled via αNIG. We refer the reader to Barndorff-Nielsen (1998) for more details on the NIG distribution.

Using maximum likelihood estimation (MLE),22We applied the built-in function “nigFit” in the R library “fBasics” for the MLE with the moment estimated parameters as initial values. we obtain the parameters for the NIG distribution fitted to the log returns of the ICAP price data as reported in Table 1. We note that these are rescaled for a NIG distribution on a daily time scale, which means that we divide the estimates for μNIG and δNIG by four. Recall that we have a trade roughly every fourth day, and we have thus supposed a sampling interval of four days for simplicity.

In Figure 2, we see the excellent fit of the NIG distribution to the ICAP log return data.

Figure 3: Time series for the daily wind production in DK1
Figure 3: Time series for the daily wind production in DK1. Data is provided by Montel, collected from the Nord Pool database.

To conduct a case study on the production rate R(t), we consider a data set of actual wind production in Denmark (the area DK1). We have available daily wind power production figures (in megawatt hours) from January 1, 2015 to November 16, 2016. The time series of the 685 production data points is plotted in Figure 3.

We recall from (2.3) that ln(R(t))=U(t), where U is an OU process. On a discrete time scale, an OU process is an AR(1) time series, and by simple linear regression we estimate the speed of mean reversion to be α^=0.51, the mean level μ^=10.3 and the standard deviation σ^u=0.71.

Figure 4: Empirical density for the residuals of the daily wind production in DK1
Figure 4: Empirical density for the residuals of the daily wind production in DK1 with the fitted normal distribution (dashed line) on a logarithmic frequency scale.

In Figure 4, we show the fitted normal distribution to the residuals of the regression analysis together with the empirical density on a logarithmic frequency scale. The normal model is far from perfect, but it still captures the distribution of the wind variation in both tails and the center reasonably well. In this analysis, we have ignored any potential seasonal effects in the production data.

With both the production and the TGC price data models at hand, we can analyze the optimal management of green certificates. The characteristic function of the NIG distribution is (see, for example, Schoutens 2003)

  ϕ(u)=exp(iuμNIG-δNIG(αNIG2-(βNIG+iu)2-αNIG2-βNIG2)).  

Thus, we let Y in (3.1) have the Lévy triplet (γ,0,ν(dξ)), where

  γ=2δNIGαNIGπ01sinh(βNIGx)K1(αNIGx)dx+μNIG  

and

  ν(dξ)=δNIGαNIGπexp(βNIGξ)K1(αNIG|ξ|)|ξ|dξ,  

and the modified Bessel function of the third kind Kv(z), with index v, is given by

  Kv(z)=120uv-1exp(-12z(u+u-1))du,z>0.  

Since the NIG distribution has finite moments, the condition in Lemma 3.1 holds and we obtain

  aX=lnϕ(-i)=μNIG-δNIG(αNIG2-(βNIG+1)2-αNIG2-βNIG2).   (4.1)

The NIG parameters in Table 1 provide us with the estimate a^X=-0.000343 and therefore aX-r<0 for any nonnegative interest rates. Note that (3.2) holds for any k<16.6 since

  k+β^NIG-12α^NIG=k-16.6,  

which is less than zero as long as k<16.6. Hence, by Theorem 3.8, the optimal strategy is to use the control A~1, and the optimal value is given by

  V(t,x,ϱ,z)=xz+cxH(t,ϱ,T).  

From the definition of H in (3.11) and Lemma 3.3(i), after a change of variables, we obtain

  V(t,x,ϱ,z)=xz+cx0T-te-ruh(u,ϱ)du,  

with (slightly abusing the notation)

  h(u,ϱ)=exp(ln(ϱ)e-αu+μα(1-e-αu)+σu22α(1-e-2αu)+aXu).  

Note that, when u becomes large,

  h(u,ϱ)exp(aXu+2μ+σu22α),  

which in our empirical example will tend to zero with an exponential rate aX<0. Therefore, the integral

  0T-texp(-ru)h(u,ϱ)du  

will tend to a constant when T-t is large. The current state, ϱ, of the wind production will, in the long term, have no effect on the current value V. This is of course not unreasonable, taking into account that our production rate is modeled as a stationary process. In fact, the half-life of the production rate is approximately 1.3 days, so it will quickly reach a stationary level.

Figure 5: The integral ∫0T-texp⁡(-r⁢u)⁢h⁢(u,ϱ)⁢d⁢u as a function of T-t.
Figure 5: The integral 0T-texp(-ru)h(u,ϱ)du as a function of T-t.

In Figure 5, we plot the integral term 0T-texp(-ru)h(u,ϱ)du as a function of T-t using the estimated parameters as input. This essentially describes the time value of the value function V. We set ln(ϱ)=μ/α, which is the stationary mean of U(t), while the annual discount rate is r=7%. After a ten-year horizon, we see that the integral has still not reached its asymptote.

In Figure 6, we zoom in on the integral value for T-t small, ie, in the range 0–10 days. This will be in the closing days of the market.

Figure 6: The integral ∫0T-texp⁡(-r⁢u)⁢h⁢(u,ϱ)⁢d⁢u as a function of T-t small.
Figure 6: The integral 0T-texp(-ru)h(u,ϱ)du as a function of T-t small.

The dotted lines illustrate the value if ln(ϱ) is 50% below (dotted line below the complete) or 50% above (dotted line above the complete) the stationary mean. There are significant effects on the time value for small T-t. Note also that the integral is close to linear in T-t, while in Figure 5 it is concave overall.

For a particular producer of renewable power qualifying for TGCs, data for the production from the power plant should of course be used in order to model the production rate. In our case study, we have used wind power from Denmark only as a proxy. The OU model can admittedly be improved. As mentioned earlier, higher-order autoregressive models (or CARMA models in continuous time) could be suitable, as wind, for example, is typically well explained by such models (see, for example, Benth and Saltyte Benth 2013). Lévy processes may also be called for when modeling residuals, as we may also observe leptokurtic behavior in such data. The main drawback with higher-order autoregressive models in our context is the optimization problem, which in a Markovian context will depend on unobserved states of the autoregressive model. This is the reason for choosing a simpler dynamics for the production rate in this paper. The TGC price dynamics seems very well suited for the purpose. However, we may also improve our analysis by performing an estimation that takes into account varying time steps between observations.

5 Proofs and intermediate results

In this section, we have collected proofs of most of the results in the paper, along with necessary intermediate results.

5.1 Proof of Proposition 3.2

The value function is defined in (2.1) as

  V(t,x,ϱ,z):=supA𝒜(t)𝔼[tTe-r(s-t)X(s)dA(s)|𝑿(t)=(x,ϱ,z)].  

By Bellman’s principle of optimality, we have, for tθT,

  0=supA𝒜(t)𝔼[tθe-rsX(s)dA(s)+e-rθV(θ,X(θ),R(θ),Z(θ))-e-rtV(t,x,ϱ,z)|𝑿(t)=(x,ϱ,z)].   (5.1)

By Ito’s formula, we obtain

  tθd(e-rsV(s,X(s),R(s),Z(s))) =tθe-rsVt(s,X(s),R(s),Z(s))-re-rsV(s,X(s),R(s),Z(s))ds  
      +tθe-rsVx(s,X(s),R(s),Z(s))aXX(s)ds  
      +tθe-rsVϱ(s,X(s),R(s),Z(s))aR(R(s))R(s)ds  
      +tθe-rsVz(s,X(s),R(s),Z(s))cR(s)ds  
      +tθe-rsVx(s,X(s),R(s),Z(s))σYX(s)dBY(s)  
      +tθe-rsVϱ(s,X(s),R(s),Z(s))σuR(s)dBu(s)  
      -tθe-rsVz(s,X(s),R(s),Z(s))dA(s)  
      +tθe-rsVxx(s,X(s),R(s),Z(s))12σY2X2(s)ds  
      +tθe-rsVϱϱ(s,X(s),R(s),Z(s))12σu2R2(s)ds  
      +tθe-rsVxϱ(s,X(s),R(s),Z(s))σYσuρX(s)R(s)ds  
      +tθe-rs[V(s,X(s-)eξ,R(s-),Z(s-))-V(s,X(s-),R(s-),Z(s-))]N~(ds,dξ)  
      +tθe-rs[V(s,X(s)+X(s-)(eξ-1),R(s),Z(s))-V(s,X(s),R(s),Z(s))-X(s)(eξ-1)]ν(dξ)ds.  

By the conditions in the proposition, we get, using Bellman’s principle,

  supA𝒜(t)𝔼[tθe-rs(X(s)-Vz)dA(s)+tθe-rs[Vt+V-rV]ds|𝑿(t)=(x,ϱ,z)]=0,   (5.2)

where the operator L is given in (3.2). Clearly, (5.2) is satisfied by the HJB equation:

  max(Vt+V-rV,-Vz+x)=0.  

5.2 Proof of Lemma 3.3

We have that

  X(s)R(s)=xeY(s)eU(s).  

Note that we can write the dynamics of Y as

  dY(v)=dY0(v)+σYdBY(v),  

where

  dY0(v):=γdv+|ξ|<1ξN~(dv,dξ)+|ξ|1ξN(dv,dξ),  

using the representation

  Y(s)=Y0(s)+σYtsdBY(v).  

Here, Y0 has the characteristic triplet (γ,0,ν(dξ)). Furthermore, due to the correlation between BY and Bu, we have

  Bu(s)=ρBY(s)+1-ρ2W(s),   (5.3)

where W(s) is another Brownian motion, independent of BY(s). It follows that

  𝔼[X(s)R(s)𝑿(t)=(x,ϱ,z)] =xexp(ln(ϱ)e-α(s-t)+μα(1-e-α(s-t)))  
      ×𝔼[exp(tsσYdBY(v)+tsσue-α(s-t)dBu(v))|𝑿(t)=(x,ϱ,z)]  
      ×𝔼[exp(Y0(s))𝑿(t)=(x,ϱ,z)].   (5.4)

For the first expectation, we have, by using (5.3),

  tsσYdBY(v)+tsσue-α(s-t)dBu(v)=ts(σue-α(s-v)ρ+σY)dBY(v)+tsσue-α(s-v)1-ρ2dW(v).   (5.5)

Both integrals in (5.5) have zero expectation, and

  Var(ts(σue-α(s-v)ρ+σY)dBY(v)) =ts(σue-α(s-v)ρ+σY)2dv  
    =σu22α(1-e-2α(s-t))+2σuσYα(1-e-α(s-t))+σY2(s-t)   (5.6)

and

  Var(tsσue-α(s-v)1-ρ2dW(v)) =tsσu2e-2α(s-v)(1-ρ2)dv  
    =σu2(1-ρ2)2α(1-e-2α(s-t)),   (5.7)

where we have used the Ito isometry. Consequently, (5.5) is normally distributed, with mean zero and variance being the sum of (5.2) and (5.2). It follows that

  𝔼[exp(σYtsdBY(v)+tsσue-α(s-v)dBu(v))|𝑿(t)=(x,ϱ,z)]=exp[σu24α(1-e-2α(s-t))(2-ρ2)+σuσYα(1-e-α(s-t))+12σY2(s-t)].  

We now turn to the second expectation in (5.2). By the Lévy–Kinchtine formula, we obtain

  𝔼[eiuY0(s)𝑿(t)=(x,ϱ,z)]=ϕY0(s-t)(u),  

where ϕY0 is the characteristic function for Y0(1). From the assumption on ϕ, which is inherited by ϕY0, it follows that, with u=-i,

  𝔼[eY0(s)𝑿(t)=(x,ϱ,z)]=ϕY0(s-t)(-i).  

Hence,

  𝔼[X(s)R(s)𝑿(t)=(x,ϱ,z)] =xexp[ln(ϱ)e-α(s-t)+μα(1-e-α(s-t))+σu24α(1-e-2α(s-t))(2-ρ2)  
            +σuσYα(1-e-α(s-t))+12σY2(s-t)+(lnϕY0(-i))(s-t)]  
    =xexp[ln(ϱ)e-α(s-t)+μα(1-e-α(s-t))+σu24α(1-e-2α(s-t))(2-ρ2)  
                    +σuσYα(1-e-α(s-t))+(lnϕ(-i))(s-t)].  

By Lemma 3.1, claim (i) follows. For (ii), we obtain

  𝔼[X(T)𝑿(t)(x,ϱ,z)] =x𝔼[e(Y0(T))𝑿(t)=(x,ϱ,z)]𝔼[exp(tTσYdBY(v))|𝑿(t)=(x,ϱ,z)]  
    =xϕ(T-t)(-i)exp(12σY2(T-t))  
    =xexp(lnϕY0(-i)+12σY2)(T-t)  
    =xexp(lnϕ(-i))(T-t)  
    =xexp(aX(T-t)).  

For (iii), we have

  𝔼[X(T)R(s)𝑿(t)=(x,ϱ,z)] =xexp(ln(ϱ)e-α(s-t)+μα(1-e-α(s-t)))  
    ×𝔼[exp(tTσYdBY(v)+tsσue-α(s-t)dBu(v))|𝑿(t)=(x,ϱ,z)]  
    ×𝔼[exp(Y0(T))𝑿(t)=(x,ϱ,z)].  

Similarly to (5.5), we get

  tTσYdBY(v)+tsσue-α(s-t)dBu(v) =sTσYdBY(v)+tsσYdBY(v)+tsσue-α(s-t)dBu(v)  
    =sTσYdBY(v)+ts(σue-α(s-v)ρ+σY)dBY(v)  
      +tsσue-α(s-v)1-ρ2dW(v).  

The first integral has expectation zero and variance σY2(T-s). Thus,

    𝔼[exp(tTσYdBY(v)+tsσue-α(s-t)dBu(v))|𝑿(t)=(x,ϱ,z)]  
        =exp[σu24α(1-e-2α(s-t))(2-ρu2)+σuσYα(1-e-α(s-t))+12σY2(T-t)],  

as in the proof of (i). Hence,

  𝔼[X(T)R(s)𝑿(t)=(x,ϱ,z)] =xexp[ln(ϱ)e-α(s-t)+μα(1-e-α(s-t))+σu24α(1-e-2α(s-t))(2-ρu2)  
             +σuσYα(1-e-α(s-t))+12σY2(T-t)lnϕ(-i)(T-t)]  
    =xexp[ln(ϱ)e-α(s-t)+μα(1-e-α(s-t))+σu24α(1-e-2α(s-t))(2-ρu2)  
             +σuσYα(1-e-α(s-t))+lnϕ(-i)(T-t)].  

The result follows.

5.3 Calculation of the derivatives of ϕ in terms of h with the control A~1

Define

  h:=h(t,ϱ,s):=1x𝔼[X(s)R(s)𝑿(t)=(x,ϱ,z)],  
  H(t,ϱ,T):=tTe-r(s-t)h(t,ϱ,s)ds.  

For the control A~1, we have

  Φ(t,x,ϱ,z)=xz+cxH(t,ϱ,T).   (5.8)

By elementary differentiation, we get

  Φx =z+ctTe-r(s-t)hds,  
  Φxx =0,  
  Φϱ =cxHϱ=cxtTe-r(s-t)hϱds,  
  Φxϱ =cHϱ=ctTe-r(s-t)hϱds,  
  Φϱϱ =cxHϱϱ=cxtTe-r(s-t)hϱϱds,  
  Φt =cxHt=cxtTe-r(s-t)(ht+rh)ds-cxh(t,ϱ,t),  
  Φz =x.  

The partial derivatives of h(t,ϱ,s) are given by

  hϱ =1ϱhe-α(s-t),  
  hϱϱ =1ϱ2h(e-2α(s-t)-e-α(s-t)),  
  ht =M(s,t)h,  

where

  M(s,t)=e-α(s-t)[αln(ϱ)-μ-σuσY]+12σu2(ρ2-2)e-2α(s-t)-aX.   (5.9)

5.4 Calculation of the derivatives of ϕ and h~ with the control A~2

Define

  h~:=h~(t,ϱ,s):=1x𝔼[X(T)R(s)𝑿(t)=(x,ϱ,z)],  
  H~(t,ϱ,T):=tTe-r(s-t)h~(t,ϱ,s)ds.  

For the control A~2, we have

  Φ(t,x,ϱ,z)=xze(aX-r)(T-t)+cxH~(t,ϱ,T).   (5.10)

By elementary differentiation, we get

  Φx =ze(aX-r)(T-t)+ctTe-r(T-t)h~ds,  
  Φxx =0,  
  Φϱ =cxH~ϱ=cxtTe-r(T-t)h~ϱds,  
  Φxϱ =cH~ϱ=ctTe-r(T-t)h~ϱds,  
  Φϱϱ =cxH~ϱϱ=cxtTe-r(s-t)h~ϱϱds,  
  Φt =xze(aX-r)(T-t)(r-aX)+cxH~t  
    =xze(aX-r)(T-t)(r-aX)  
      +cxtTe-r(T-t)(h~t+rh~)ds-cxϱe(aX-r)(T-t),  
  Φz =xe(aX-r)(T-t).  

The partial derivatives of h~(t,ϱ,s) are given by

  h~ϱ =1ϱh~e-α(s-t),  
  h~ϱϱ =1ϱ2h~(e-2α(s-t)-e-α(s-t)),  
  h~t =M(s,t)h~,  

where

  M(s,t)=e-α(s-t)[αln(ϱ)-μ-σuσY]+12σu2(ρ2-2)e-2α(s-t)-aX.  

5.5 Proof of Lemma 3.7

We start with a result stated in, for example, Cont and Tankov (2004). Let ψin for i=1,2 be any simple predictable function. Then, the process

  θtθ{0}ψin(s,ξ)dN~(ds,dξ)   (5.11)

is a square integrable martingale that verifies the isometry formula.

Define

  Φ1(s,X(s),R(s),Z(s)) :=X(s)Z(s)+cX(s)H(s,R(s),T)   (5.12)
and
  Φ2(s,X(s),R(s),Z(s)) :=X(s)Z(s)e(aX(2)-r)(T-s)+cX(s)H~(s,R(s),T),   (5.13)

where aX(i) is associated with Φi. Define

  ψi(s,ξ) :=e-rs[Φi(s,X(s)eξ,R(s),Z(s))-Φi(s,X(s),R(s),Z(s))]  
    =e-rs(eξ-1)Φi(s,X(s),R(s),Z(s)).  

If

  𝔼[tθ{0}|ψi(s,ξ)|2ν(dξ)ds|𝑿(t)=(x,ϱ,z)]<   (5.14)

holds for i=1,2, then there exists a sequence (ψin) of simple predictable functions such that (5.11) converges, in L2(𝑷), to a process

  θtθ{0}ψi(s,ξ)dN~(ds,dξ).   (5.15)

The limiting process (5.15) is also a square integrable martingale that verifies the isometry formula (see, for example, Cont and Tankov 2004; Ikeda and Watanabe 1981). Thus, the martingale property of (3.19) follows if (5.14) holds. By Fubini’s Theorem, and since r>0, we have

  𝔼[tθ{0}|ψi(s,ξ)|2ν(dξ)ds|𝑿(t)=(x,ϱ,z)]{0}|eξ-1|ν(dξ)𝔼[tθ|Φi(s,X(s),R(s),Z(s))|2ds|𝑿(t)=(x,ϱ,z)].   (5.16)

The first integral is finite by condition (3.2). For notational convenience, define

  D[μ,α,σu,σY,ρ](τ-s):=μα(1-e-α(τ-s))+σu24α(1-e-2α(τ-s))(2-ρ2)+σuσYα(1-e-α(τ-s)).  

Then,

  Φ1(s,X(s),R(s),Z(s)) :=X(s)Z(s)+csTX(s)R(s)exp[-α(τ-s)]exp((aX(1)-r)(τ-s))D[μ,α,σu,σY,ρ](τ-s)dτ,   (5.17)

and

  Φ2(s,X(s),R(s),Z(s)) :=exp((aX(2)-r)(T-s))[X(s)Z(s)+csTX(s)R(s)exp[-α(τ-s)]D[μ,α,σu,σY,ρ](τ-s)dτ].   (5.18)

Recall that, by assumption, aX(1)-r<0 and aX(2)-r>0. Hence,

  |Φ1(s,X(s),R(s),Z(s))| X(s)Z(s)+csTX(s)R(s)exp[-α(τ-s)]D[μ,α,σu,σY,ρ](τ-s)dτ  
    |Φ2(s,X(s),R(s),Z(s))|.   (5.19)

It follows from (5.16) that (5.14) holds for i=1,2 if

  𝔼[tθ|Φ2(s,X(s),R(s),Z(s))|2ds|𝑿(t)=(x,ϱ,z)]   (5.20)

is finite. We obtain

  |Φ2(s,X(s),R(s),Z(s))|2 =e2(aX(2)-r)(T-s)[X(s)Z(s)+csTX(s)R(s)exp[-α(τ-s)]D[μ,α,σu,σY,ρ](τ-s)dτ]2  
    =e2(aX(2)-r)(T-s)[X2(s)Z2(s)+2csTX2(s)Z(s)R(s)exp[-α(τ-s)]D[μ,α,σu,σY,ρ](τ-s)dτ  
            +c2sTX(s)R(s)exp[-α(τ1-s)]D[μ,α,σu,σY,ρ](τ1-s)dτ1sTX(s)R(s)exp[-α(τ2-s)]D[μ,α,σu,σY,ρ](τ2-s)dτ2]  
    =e2(aX(2)-r)(T-s)[X2(s)Z2(s)+2csTX2(s)Z(s)R(s)exp[-α(τ-s)]𝟏(R(s)1)D[μ,α,σu,σY,ρ](τ-s)dτ  
            +2csTX2(s)Z(s)R(s)exp[-α(τ-s)]𝟏(R(s)<1)D[μ,α,σu,σY,ρ](τ-s)dτ  
            +c2sTsT[X2(s)R(s)e-α(τ1+τ2-2s)𝟏(R(s)1)D[μ,α,σu,σY,ρ](τ1-s)D[μ,α,σu,σY,ρ](τ2-s)]dτ1dτ2  
            +c2sTsT[X2(s)R(s)e-α(τ1+τ2-2s)𝟏(R(s)<1)D[μ,α,σu,σY,ρ](τ1-s)D[μ,α,σu,σY,ρ](τ2-s)]dτ1dτ2].  

Taking conditional expectation, using Fubini’s Theorem and that e-α()(0,1), we obtain

    𝔼[|Φ2(s,X(s),R(s),Z(s))|2𝑿(t)=(x,ϱ,z)]  
    e2(aX(2)-r)(T-s)[𝔼[X2(s)Z2(s)𝑿(t)=(x,ϱ,z)]  
            +2csT𝔼[X2(s)Z(s)R(s)𝑿(t)=(x,ϱ,z)]D[μ,α,σu,σY,ρ](τ-s)dτ  
            +2csT𝔼[X2(s)Z(s)𝑿(t)=(x,ϱ,z)]D[μ,α,σu,σY,ρ](τ-s)dτ  
            +c2sTsT(𝔼[X2(s)R(s)𝑿(t)=(x,ϱ,z)]D[μ,α,σu,σY,ρ](τ1-s)D[μ,α,σu,σY,ρ](τ2-s))dτ1dτ2  
            +c2sTsT(𝔼[X2(s)𝑿(t)=(x,ϱ,z)]D[μ,α,σu,σY,ρ](τ1-s)D[μ,α,σu,σY,ρ](τ2-s))dτ1dτ2].  

Note that

  D[μ,α,σu,σY,ρ](τ-s){exp(μα+σu24α+σuσYα)if μ0,exp(σu24α+σuσYα)if μ<0.  

Define

  D¯[μ,α,σu,σY,ρ]:=max[exp(μα+σu24α+σuσYα),exp(σu24α+σuσYα)].  

Then,

  D[μ,α,σu,σY,ρ](τ-s)D¯[μ,α,σu,σY,ρ].  

We obtain

  𝔼[|Φ2(s,X(s),R(s),Z(s))|2𝑿(t)=(x,ϱ,z)] e2(aX(2)-r)(T-s)[𝔼[X2(s)Z2(s)𝑿(t)=(x,ϱ,z)]  
            +2c(T-s)D¯[μ,α,σu,σY,ρ]𝔼[X2(s)Z(s)R(s)𝑿(t)=(x,ϱ,z)]  
            +2c(T-s)D¯[μ,α,σu,σY,ρ]𝔼[X2(s)Z(s)𝑿(t)=(x,ϱ,z)]  
            +c2(T-s)2D¯2[μ,α,σu,σY,ρ]𝔼[X2(s)R(s)𝑿(t)=(x,ϱ,z)]  
                    +c2(T-s)2D¯2[μ,α,σu,σY,ρ]𝔼[X2(s)𝑿(t)=(x,ϱ,z)]].   (5.21)

Recall that

  Z(s)=cP(s)-A(s),where P(s)=0sR(v)dv.  

Since A is nondecreasing and A(s)cP(s), we have

  Z(s)ctsR(v)dv+Z(t),Z2(s)[ctsR(v)dv+Z(t)]2=c2tsR(u)Z(u)du+2cZ(t)tsR(v)dv+Z2(t).}   (5.22)

We now calculate the conditional expectations in (5.5). In the calculations below, we will use

  𝔼[𝑿(t)]=𝔼[t].  

For 𝔼[X2(s)𝑿(t)=(x,ϱ,z)]

Similarly to the proof of Lemma 3.3(ii), we get

  𝔼[X2(s)t] =X2(t)𝔼[e2Y0(s)t]𝔼[exp(ts2σYdBY(v))|t]  
    =X2(t)exp(ln(ϕ(-2i))(s-t)).   (5.23)

Since X(t)=x, we have

  𝔼[X2(s)𝑿(t)=(x,ϱ,z)]=x2exp(ln(ϕ(-2i))(s-t)).   (5.24)

This is clearly positive and finite for s[t,T].

For 𝔼[X2(s)R(s)𝑿(t)=(x,ϱ,z)]

Similarly to the proof of Lemma 3.3(i), we obtain

  𝔼[X2(s)R(s)t] =𝔼[X2(t)e2Y(s)eU(s)t]  
    =X2(t)R(t)exp[-α(s-t)]exp(μα(1-e-α(s-t)))  
      ×𝔼[exp(ts2σYdBY(v)+tsσue-α(s-t)dBu(v))|t]𝔼[e2Y0(s)t]  
    =X2(t)R(t)exp[-α(s-t)]D[μ,α,σu,2σY,ρ](s-t)exp[ln(ϕ(-2i))(s-t)]  
    =X2(t)R(t)exp[-α(s-t)]𝟏(R(t)1)×D[μ,α,σu,2σY,ρ](s-t)exp[ln(ϕ(-2i))(s-t)]  
      +X2(t)R(t)exp[-α(s-t)]𝟏(R(t)<1)D[μ,α,σu,2σY,ρ](s-t)exp[ln(ϕ(-2i))(s-t)]  
    X2(t)(R(t)+1)D¯[μ,α,σu,2σY,ρ]exp[ln(ϕ(-2i))(s-t)].   (5.25)

Since X(t)=x,R(t)=ϱ, we have

  𝔼[X2(s)R(s)𝑿(t)=(x,ϱ,z)]x2(ϱ+1)D¯[μ,α,σu,2σY,ρ]exp[ln(ϕ(-2i))(s-t)].   (5.26)

This is clearly positive and finite for s[t,T].

For 𝔼[X2(s)Z(s)𝑿(t)=(x,ϱ,z)]

By (5.22), we obtain

  𝔼[X2(s)Z(s)t]𝔼[X2(s)ctsR(v)dv|t]+Z(t)𝔼[X2(s)t].   (5.27)

The last term is positive and finite by (5.24). Thus, by Fubini’s Theorem and the tower property, we obtain

  𝔼[X2(s)Z(s)t] cts𝔼[X2(s)R(v)t]dv  
    =cts𝔼[R(v)𝔼[X2(s)v]t]dv  
    =cts𝔼[X2(v)R(v)t]exp(ln(ϕ(-2i))(s-v))dv,   (5.28)

where we have used (5.5) with t replaced by v. By (5.5), we have

  𝔼[X2(v)R(v)t]X2(t)(R(t)+1)D¯[μ,α,σu,2σY,ρ]exp[ln(ϕ(-2i))(v-t)].   (5.29)

Hence,

  𝔼[X2(s)Z(s)t] ctsX2(t)(R(t)+1)D¯[μ,α,σu,2σY,ρ]exp(ln(ϕ(-2i))(s-t))dv  
    =cX2(t)(R(t)+1)(s-t)D¯[μ,α,σu,2σY,ρ]exp(ln(ϕ(-2i))(s-t)).   (5.30)

Since X(t)=x and R(t)=ϱ, we have

  𝔼[X2(s)Z(s)𝑿(t)=(x,ϱ,z)]cx2(ϱ+1)(s-t)D¯[μ,α,σu,2σY,ρ]exp(ln(ϕ(-2i))(s-t)).   (5.31)

This is clearly positive and finite for s[t,T].

For 𝔼[X2(s)Z(s)R(s)𝑿(t)=(x,ϱ,z)]

By (5.22), (5.26), Fubini’s Theorem and the tower property, we have

  𝔼[X2(s)Z(s)R(s)t] 𝔼[X2(s)R(s)tsR(v)dv|t]+Z(t)𝔼[X2(s)R(s)t]  
    ts𝔼[R(v)𝔼[X2(s)R(s)v]t]dv.   (5.32)

As in (5.5), we have

  𝔼[X2(s)R(s)v]=X2(v)R(v)exp[-α(s-v)]D[μ,α,σu,2σY,ρ](s-v)exp[ln(ϕ(-2i))(s-v)].  

Hence,

  𝔼[X2(s)Z(s)R(s)t] ts𝔼[X2(v)R(v)(1+e-α(s-v))t]D[μ,α,σu,2σY,ρ](s-v)exp[ln(ϕ(-2i))(s-v)]dv.  
    =ts(𝔼[X2(v)R(v)(1+e-α(s-v))𝟏(R(v)1)t]D[μ,α,σu,2σY,ρ](s-v)exp[ln(ϕ(-2i))(s-v)])dv.  
      +ts(𝔼[X2(v)R(v)(1+e-α(s-v))𝟏(R(v)<1)t]D[μ,α,σu,2σY,ρ](s-v)exp[ln(ϕ(-2i))(s-v)])dv.  
    ts𝔼[X2(v)R(v)2t]D¯[μ,α,σu,2σY,ρ]exp[ln(ϕ(-2i))(s-v)]dv  
      +ts𝔼[X2(v)t]D¯[μ,α,σu,2σY,ρ]exp[ln(ϕ(-2i))(s-v)]dv.  

The last term is positive and finite by (5.24). For 𝔼[X2(v)R(v)2t], we obtain, similarly to the derivation of (5.5),

  𝔼[X2(v)R(v)2t] =X2(t)R(t)2exp[-α(v-t)]D[2μ,α,2σu,2σY,ρ](v-t)exp[ln(ϕ(-2i))(v-t)]  
    X2(t)R2(t)𝟏(R(v)1)D[2μ,α,2σu,2σY,ρ](v-t)exp[ln(ϕ(-2i))(v-t)]  
      +X2(t)𝟏(R(v)<1)D[2μ,α,2σu,2σY,ρ](v-t)exp[ln(ϕ(-2i))(v-t)]  
    X2(t)(R2(t)+1)D¯[2μ,α,2σu,2σY,ρ]exp[ln(ϕ(-2i))(v-t)].  

Hence,

  𝔼[X2(s)Z(s)R(s)t] tsX2(t)(R2(t)+1)D¯[2μ,α,2σu,2σY,ρ]  
      ×D¯[μ,α,σu,2σY,ρ]exp[ln(ϕ(-2i))(s-t)]  
    =X2(t)(R2(t)+1)(s-t)D¯[2μ,α,2σu,2σY,ρ]  
      ×D¯[μ,α,σu,2σY,ρ]exp[ln(ϕ(-2i))(s-t)].  

Since X(t)=x,R(t)=ϱ, we have

  𝔼[X2(s)Z(s)R(s)𝑿(t)=(x,ϱ,z)]x2(t)(ϱ2(t)+1)(s-t)D¯[2μ,α,2σu,2σY,ρ]D¯[μ,α,σu,2σY,ρ]exp[ln(ϕ(-2i))(s-t)].   (5.33)

This is clearly positive and finite for s[t,T].

For 𝔼[X2(s)Z2(s)t]

Again, by (5.22) and Fubini’s Theorem, we have

  𝔼[X2(s)Z2(s)t] c2ts𝔼[X2(s)R(v)Z(v)t]dv  
      +2cZ(t)ts𝔼[X2(s)R(v)t]dv  
      +Z2(t)ts𝔼[X2(s)t]dv.  

The last two terms are found to be positive and finite via the tower property and (5.26) and (5.24). Hence,

  𝔼[X2(s)Z2(s)t] c2ts𝔼[X2(s)R(v)Z(v)t]dv  
    =c2ts𝔼[R(v)Z(v)𝔼[X2v]t]dv  
    =c2tsexp[ln(ϕ(-2i))(s-v)]𝔼[X2(v)R(v)Z(v)t]dv  
    c2X2(t)(R2(t)+1)D¯[2μ,α,2σu,2σY,ρ]D¯[μ,α,σu,2σY,ρ]exp[ln(ϕ(-2i))(s-t)]ts(v-t)dv  
    =c4X2(t)(R2(t)+1)(s-t)2D¯[2μ,α,2σu,2σY,ρ]D¯[μ,α,σu,2σY,ρ]exp[ln(ϕ(-2i))(s-t)].  

We get

  𝔼[X2(s)Z2(s)𝑿(t)=(x,ϱ,z)] c4x2(ϱ2+1)(s-t)2D¯[2μ,α,2σu,2σY,ρ]D¯[μ,α,σu,2σY,ρ]exp[ln(ϕ(-2i))(s-t)].   (5.34)

This is clearly positive and finite for s[t,T]. All the terms in (5.5) are finite and continuous in s. The square integrability of (5.16) follows, and (5.14) is indeed finite. It follows that the process

  θtθe-rs[Φ(s,X(s)eξ,R(s),Z(s))-Φ(s,X(s),R(s),Z(s))]N~(ds,dξ)  

is a martingale. To see that the processes

  θ tθe-rsΦx(s,X(s),R(s),Z(s))σYX(s)dBY(s),  
  θ tθe-rsΦϱ(s,X(s),R(s),Z(s))σuR(s)dBu(s)  

are martingales, note that

  tθe-rsΦix(s,X(s),R(s),Z(s))σYX(s)dBY(s)=tθe-rsΦi(s,X(s),R(s),Z(s))σYX(s)dBY(s)  

and

  R(s)Φiϱ(s,X(s),R(s),Z(s))Φi(s,X(s),R(s),Z(s)).   (5.35)

Hence, the martingale property follows, as for (3.19), directly from the finiteness of (5.20). The statement of the lemma follows.

6 Conclusions

We provided a valuation model for the income of selling TGCs, formulated as a singular stochastic control problem. Our model takes into account the production rate of renewable energy from a “typical” plant, the market price of TGCs and the cumulative number of certificates sold. We assumed that the production rate dynamics is given by an exponential OU process and that the logarithm of the TGC price is a Lévy process.

We found an explicit solution to the control problem, including the optimal selling strategy. The optimal value is easily calculated by Theorem 3.8 via Lemma 3.3 as soon as the characteristic function of the logarithmic price is known.

Furthermore, we conducted an empirical analysis on ICAP TGC prices, where a NIG-distributed Lévy process appeared to be a highly appropriate choice to describe the dynamics. We also fitted the production rate model to a sample of realized wind power production in Denmark. The model functioned reasonably well in explaining the dynamics. Given the fitted models, we observed that the optimal strategy for a wind power producer in this market is to immediately sell the TGCs granted. We analyzed further the time value of TGC for a producer and how the mean reversion in renewable production affects the value of TGCs to be sold.

The models and analysis in this paper can be extended in several ways. First, it may be highly unrealistic to assume a constant discount rate r over a decade. An extension could be to assume a stochastic (spot) discount rate of Markovian type (for example, the Vasicek model). This would create a slightly more involved optimization problem. We could also investigate the price process model for TGCs, as government interventions in the TGC market that regulate demand for certificates or structural changes in the power market (eg, connecting cables, construction of new power plants) may affect the TGC price dynamics in various ways. A possible extension here would be to introduce multifactor models, which yield more flexibility in capturing different stochastic variations in mean, drift and volatility. This would lead to an even more demanding stochastic control problem to analyze.

Declaration of interest

We are grateful to Montel AS for providing us with data. F. E. Benth acknowledges financial support from the project FINEWSTOCH, funded by the Norwegian Research Council. The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.

Acknowledgements

We thank the anonymous referee for criticism and suggestions significantly improving the presentation of the paper.

References

Amundsen, E. S., and Mortensen, J. B. (2001). The Danish green certificate system: some simple analytical results. Energy Economics 23(6), 489–509 (http://doi.org/fdjz6j).

Aune, F. R., Dalen, H. M., and Hagem, C. (2012). Implementing the EU renewable target through green certificate markets. Energy Economics 34(4), 992–1000 (http://doi.org/b3tc8k).

Barndorff-Nielsen, O. E. (1998). Processes of normal inverse Gaussian type. Finance & Stochastics 2(1), 41–68 (http://doi.org/c4jr7k).

Benth, F. E., and Saltyte Benth, J. (2013). Modeling and Pricing in Financial Markets for Weather Derivatives. Advanced Series on Statistical Science and Applied Probability, Volume 17. World Scientific, Singapore.

Benth, F. E., Saltyte Benth, J., and Koekebakker, S. (2008). Stochastic Modelling of Electricity and Related Markets. Advanced Series on Statistical Science and Applied Probability, Volume 11. World Scientific, Singapore.

Cont, R., and Tankov, P. (2004). Financial Modelling with Jump Processes. Chapman & Hall/CRC, Boca Raton, FL.

Coulon, M., Khazaei, J., and Powell, B. W. (2015). SMART-SREC: a stochastic model of the New Jersey solar renewable energy certificate market. Journal of Environmental Economics and Management 73, 13–31 (http://doi.org/b4q3).

Eriksson, A., Ghysels, E., and Wang, F. (2009). The normal inverse Gaussian distribution and the pricing of derivatives. Journal of Derivatives 16(3), 23–37 (http://doi.org/c74zg5).

Fagiani, R., and Hakvoort, R. (2014). The role of regulatory uncertainty in certificate markets: a case study of the Swedish/Norwegian market. Energy Policy 65, 608–618 (http://doi.org/b4q4).

Frogner, J. S., and Hustveit, M. (2015). An analysis of the Swedish/Norwegian electricity certificate market. Master’s Thesis, June, NTNU Norway (http://bit.ly/2ndc3d6).

Goldstein, H. S. (2010). A Green certificate market in Norway. Term Paper, Spring, Energy Economics and Policy, ETH Zurich (http://bit.ly/2mwMXKo).

Ikeda, N., and Watanabe, S. (1981). Stochastic Differential Equations and Diffusion Processes. Kodansha, Tokyo.

NVE (2013). The Swedish–Norwegian electricity certificate market annual report 2012. Norges Vassdrags- og Energidirektorat [Norwegian Water Resources and Energy Directorate]. URL: http://bit.ly/2nuWKAj.

Schoutens, W. (2003). Lévy Processes in Finance: Pricing Financial Derivatives. Wiley Series in Probability and Statistics. Wiley (http://doi.org/bkj7xq).

Unger, T., and Ahlgren, E. O. (2005). Impacts of a common green certificate market on electricity and CO2-emission markets in the Nordic countries. Energy Policy 33(16), 2152–2163 (http://doi.org/cdrd32).

Wolfgang, O., Jaehnert, S., and Mo, B. (2015). Methodology for forecasting in the Swedish–Norwegian market for el-certificates. Energy 88, 322–333 (http://doi.org/f7r3fv).

You need to sign in to use this feature. If you don’t have a Risk.net account, please register for a trial.

Sign in
You are currently on corporate access.

To use this feature you will need an individual account. If you have one already please sign in.

Sign in.

Alternatively you can request an indvidual account here: