Empirical performance of reduced-form models for emission permit prices

Abstract

The value of emission permits in environmental markets derives from the particular design features of the underlying cap-and-trade system. In this paper, we evaluate a model framework for the price dynamics of emission permits which accounts for these features in a reduced-form way. Based on permit futures and option data from the European Union Emissions Trading System, the world’s largest environmental market, we show that model variants which represent the design of the system most accurately provide the best fit to historical futures prices and achieve the best option pricing performance. Our results suggest that the specific design of cap-and-trade systems directly translates to the dynamic properties of emission permit prices, and that models tailored to environmental markets are the best choice for related pricing and risk management decisions.

Appendices

Appendix A: Dynamics of risk-neutral shortage probabilities

We derive the dynamics of the risk-neutral shortage probabilities $A_k$ in general and show that all emissions processes (3) satisfying a narrow-sense linear SDE after transformation by a strictly increasing function lead to the same class of dynamics given by (6). For notational ease we omit the k indicating the compliance period in this section.

Let $D_{t}(x_{t,T};.)$ be the cumulative density function of $x_{T}$ given $x_{t,T}$ at time t. By definition, it follows

$$\begin{aligned} A_{t}=1-D_{t}(x_{t,T};q)=:G_{t}(x_{t,T}). \end{aligned}$$

(27)

Applying Itô’s Lemma directly yields

$$\begin{aligned} { dA}_{t}=g_{t}(x_{t,T})\sigma (t,x_{t,T}) dW_t =g_{t}(G^{-1}_{t}(A_{t}))\sigma (t,G^{-1}_{t}(A_{t})) dW_t, \end{aligned}$$

(28)

where $g_{t}$ is the first derivative of $G_{t}$. Equation (28) describes the general dynamics of the risk-neutral shortage probabilities.

Aiming at simple dynamics of $A_t$, we assume that $x_{t,T}$ can be transformed by a strictly increasing function $\upsilon$ in such a way that $\upsilon (x_{t,T})$ follows a linear SDE in the narrow sense, i.e.,

$$\begin{aligned} d\upsilon (x_{t,T})=(a_1(t)\upsilon (x_{t,T})+a_2(t))dt+b(t)dW_t, \end{aligned}$$

(29)

where $a_1:(0,T)\rightarrow \mathbb {R}$ and $a_2:(0,T)\rightarrow \mathbb {R}$ are continuous, bounded functions and $b:(0,T)\rightarrow \mathbb {R}^+$ is continuous and square-integrable. This especially covers a geometric Brownian motion as proposed by Carmona and Hinz (2011) or the case that $\ln {x_{t,T}}$ follows an Ornstein-Uhlenbeck process $d\ln {x_{t,T}}=\kappa (\mu (t)-\ln {x_{t,T}}) dt+\sigma (t) dW_t$. Given $x_{t,T}$, we can write $\upsilon (x_{T})$ explicitly as

$$\begin{aligned} \upsilon (x_{T})=e^{\int _t^T a_1(s)ds}\upsilon (x_{t,T})+\int _t^T e^{\int _s^T a_1(u)du}a_2(s)ds +\int _t^T e^{\int _s^T a_1(u)du}b(s)dW_s,\nonumber \\ \end{aligned}$$

(30)

see Karatzas and Shreve (1991), pp. 360–361. Particularly, $\upsilon (x_{T})$ is normally distributed with mean

$$\begin{aligned} \mu _{\upsilon t,T}(x_{t,T})=e^{\int _t^T a_1(s)ds}\upsilon (x_{t,T})+\int _t^T e^{\int _s^T a_1(u)du}a_2(s)ds \end{aligned}$$

(31)

and standard deviation

$$\begin{aligned} \sigma _{\upsilon t,T}=\sqrt{\int _t^T e^{2\int _s^T a_1(u)du}b^2(s)ds}. \end{aligned}$$

(32)

It follows that the cumulative density function is given by

$$\begin{aligned} D_{t,T}(x_{t,T};y)=\Phi \left( \frac{\upsilon (y)-\mu _{\upsilon t,T}(x_{t,T})}{\sigma _{\upsilon t,T}}\right) , \end{aligned}$$

(33)

and we obtain

$$\begin{aligned} g_{t,T}(x_{t,T})=\Phi '\left( \frac{\upsilon (q)-\mu _{\upsilon t,T}(x_{t,T})}{\sigma _{\upsilon t,T}}\right) \frac{\mu _{\upsilon t,T}'(x_{t,T})}{\sigma _{\upsilon t,T}}. \end{aligned}$$

(34)

Inserting into (28) yields

$$\begin{aligned} \begin{aligned} { dA}_{t}&= \Phi '\left( \Phi ^{-1}(A_{t})\right) \frac{\mu _{\upsilon t,T}'(x_{t,T})\sigma (t,x_{t,T})}{\sigma _{\upsilon t,T}} dW_t \\&=\Phi '\left( \Phi ^{-1}(A_{t})\right) \frac{e^{\int _t^T a_1(s)ds}\upsilon '(x_{t,T})\sigma (t,x_{t,T})}{\sqrt{\int _t^T e^{2\int _s^T a_1(u)du}b^2(s)ds}} dW_t \end{aligned} \end{aligned}$$

(35)

for the risk-neutral shortage probability. Now observe that the dynamics of $\upsilon (x_{t,T})$ is also given by

$$\begin{aligned} d\upsilon (x_{t,T})=(\upsilon '(x_{t,T})\mu (t,x_{t,T})+\frac{1}{2}\upsilon ''(x_{t,T})\sigma ^2(t,x_{t,T}))dt +\upsilon '(x_{t,T})\sigma (t,x_{t,T})dW_t,\nonumber \\ \end{aligned}$$

(36)

applying Itô’s Lemma to (3). Comparing this to (29) yields

$$\begin{aligned} \upsilon '(x_{t,T})\sigma (t,x_{t,T})=b(t), \end{aligned}$$

(37)

and we can write (35) as

$$\begin{aligned} \begin{aligned} { dA}_{t}&=\Phi '\left( \Phi ^{-1}(A_{t})\right) \frac{e^{\int _t^T a_1(s)ds}b(t)}{\sqrt{\int _t^T e^{2\int _s^T a_1(u)du}b^2(s)ds}} dW_t \\&=\Phi '\left( \Phi ^{-1}(A_{t})\right) \frac{c(t)}{\sqrt{\int _t^T c^2(s)ds}} dW_t \end{aligned} \end{aligned}$$

(38)

with $c(t)=e^{\int _t^T a_1(s)ds}b(t)$. Carmona and Hinz (2011) show that for every continuous function $z:(0,T)\rightarrow \mathbb {R}^+$ satisfying

$$\begin{aligned} \lim _{t\rightarrow T} \int _0^t z(s) ds=\infty , \end{aligned}$$

(39)

there exists a square-integrable continuous function $c:(0,T)\rightarrow \mathbb {R}^+$ fulfilling $\frac{c^2(t)}{\int _t^T c^2(s)ds}=z(t)$. Therefore we can construct every continuous function z satisfying (39) by the choice $a_1(t)=0$ and $b(t)=c(t)$. The other way round, $c(t)=e^{\int _t^T a_1(s)ds}b(t)$ is a square-integrable continuous and positive function by the properties of $a_1$ and b, and it follows that $z(t)=\frac{c^2(t)}{\int _t^T c^2(s)ds}$ is positive and continuous and satisfies (39), see Carmona and Hinz (2011).

Altogether, for all dynamics of $x_{t,T}$ for which a strictly increasing function $\upsilon$ exists such that $\upsilon (x_{t,T})$ follows a narrow-sense linear SDE (29), the class of possible dynamics for $A_t$ is completely characterized by (6), with a continuous function $z:(0,T)\rightarrow \mathbb {R}^+$ satisfying (39). Particularly, all possible dynamics for $A_t$ in this class can be obtained by choosing an arithmetic Brownian motion (4) for the expected cumulative emissions process $x_{t,T}$.

Appendix B: Evaluation of option pricing formulas

To calculate theoretical option prices within our reduced-form model framework, Eq. (14) requires the numerical evaluation of an $m+1$-dimensional integral. While this is straightforward in the one-dimensional case ($m=0$), the following transformations may help to improve computational efficiency for models with more price components. We demonstrate these transformations for the case of $m=1$ and an additional component $R_t$ as well as for $m=2$ and $R_t=0$.

In both cases, we decompose the bivariate normal distribution as proposed by Carmona and Hinz (2011) by defining

$$\begin{aligned} \mu ^{2,c}(x_1)=\mu ^2+\frac{\nu ^{1,2}}{\nu ^{1,1}}\left( x_1-\mu ^1\right) \quad \text {and}\quad \nu ^{2,2,c}=\nu ^{2,2}-\frac{(\nu ^{1,2})^2}{\nu ^{1,1}}, \end{aligned}$$

(40)

and we can write the integral in (14) as

$$\begin{aligned}&\int _{\mathbb {R}}\int _{\mathbb {R}} \left( e^{-r(T_1-{\overline{t}})}\Phi (x_1)p_1+e^{-r(T_1-{\overline{t}})}e^{x_{2}}-K\right) ^+\nonumber \\&\quad \cdot \,\varphi \left( \mu ^{2,c}(x_1),\nu ^{2,2,c};x_2\right) dx_2 \varphi \left( \mu ^1,\nu ^{1,1};x_1\right) dx_1 \end{aligned}$$

(41)

for the case $m=1$ with an additional component $R_t$ or as

$$\begin{aligned}&\int _{\mathbb {R}}\int _{\mathbb {R}} \left( e^{-r(T_1-{\overline{t}})}\Phi (x_1)p_1+e^{-r(T_2-{\overline{t}})}\Phi (x_2)p_2-K\right) ^+ \nonumber \\&\quad \cdot \,\varphi \left( \mu ^{2,c}(x_1),\nu ^{2,2,c};x_2\right) dx_2 \varphi \left( \mu ^1,\nu ^{1,1};x_1\right) dx_1 \end{aligned}$$

(42)

for $m=2$ and $R_t=0$.

Defining $K^*(x_1)=K-e^{-r(T_1-{\overline{t}})}\Phi (x_1)p_1$, we note that the inner integral in (41) can completely be settled analytically according to

$$\begin{aligned} \begin{aligned}&\int _{\mathbb {R}}\left( e^{-r(T_1-{\overline{t}})}e^{x_{2}}-K^*(x_1)\right) ^+ \cdot \varphi \left( \mu ^{2,c}(x_1),\nu ^{2,2,c};x_2\right) dx_2 \\&=\left\{ \begin{array}{ll} e^{-r(T_1-{\overline{t}})} \int _{\mathbb {R}} e^{x_{2}} \cdot \varphi \left( \mu ^{2,c}(x_1),\nu ^{2,2,c};x_2\right) dx_2 -K^*(x_1), &{} \hbox {if }\quad K^*(x_1)\le 0; \\ \int _{\mathbb {R}}\left( e^{-r(T_1-{\overline{t}})}e^{x_{2}}-K^*(x_1)\right) ^+ \cdot \varphi (\mu ^{2,c}(x_1),\nu ^{2,2,c};x_2) dx_2, &{}\quad \hbox {otherwise.} \end{array} \right. \\&=\left\{ \begin{array}{ll} e^{-r(T_1-{\overline{t}})}\cdot e^{\mu ^{2,c}(x_1)+\frac{\nu ^{2,2,c}}{2}}-K^*(x_1), &{} \hbox {if }\quad K^*(x_1)\le 0; \\ e^{-r(T_1-{\overline{t}})}\cdot e^{\mu ^{2,c}(x_1)+\frac{\nu ^{2,2,c}}{2}}\Phi (d_1)-K^*(x_1)\Phi (d_2), &{}\quad \hbox {otherwise.} \end{array} \right. \end{aligned} \end{aligned}$$

(43)

with

$$\begin{aligned} d_1=\frac{\mu ^{2,c}(x_1)-\ln (\frac{K^*(x_1)}{e^{-r(T_1-{\overline{t}})}})}{\sqrt{\nu ^{2,2,c}}}+\sqrt{\nu ^{2,2,c}} \quad \text {and} \quad d_2=d_1-\sqrt{\nu ^{2,2,c}}. \end{aligned}$$

(44)

For (42), Carmona and Hinz (2011) show that numerical integration of the inner integral is only necessary if $0<K^*(x_1)<e^{-r(T_2-{\overline{t}})}p_2$ since outside this interval we have

$$\begin{aligned} \begin{aligned}&\int _{\mathbb {R}}\left( e^{-r(T_2-{\overline{t}})}\Phi (x_{2})p_2-K^*(x_1)\right) ^+ \cdot \varphi \left( \mu ^{2,c}(x_1),\nu ^{2,2,c};x_2\right) dx_2 \\&=\left\{ \begin{array}{ll} 0 &{}\quad \hbox {if }\quad K^*(x_1)\ge e^{-r(T_2-{\overline{t}})}p_2; \\ e^{-r(T_2-{\overline{t}})}p_2\Phi \left( \frac{\mu ^{2,c}(x_1)}{\sqrt{1+\nu ^{2,2,c}}}\right) -K^*(x_1) &{}\quad \hbox {if }\quad K^*(x_1)\le 0. \end{array} \right. \end{aligned} \end{aligned}$$

(45)

It is straightforward again to evaluate the outer integral of (41) and (42) once the inner integral is calculated.