Carbon Emission Assignment

Introduction

Energy remains an essential pillar for human livelihood for its necessity in sustainability and development of any current civilization. One of the main feature distinguishing it from most of the other commodities is the mean-reverting characteristic together with evident spikes and high volatility as seen for electricity prices. Here, we’ll look at carbon emission in Kenya. A number of works have employed the Ornstein–Uhlenbeck process to model directly the dynamics over time of different commodity spot prices under reduced-form one-factor models (Ribeiro and Hodges, 2004).

Literature Review

The Ornstein-Uhlenbeck process (Uhlenbeck, 1930) is therefore mainly employed in energy analytics to account for the mean-reverting feature. Generally, under a filtered probability space \((\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \geq 0}, \mathbb{P})\), a stochastic process \({X(t); t \geq 0}\) is said to follow a O-U process if it satisfies the SDE; \(dX_t = \theta(\mu-X_t) dt + \sigma dW_t)\) where $dX_t $ is the increment of the process \(X\) between times \(t\) and \(dt\), and \(\sigma > 0\) is the instantaneous diffusion term used as a volatility measure, \(\mu\) is the long term mean or expected value of the \(X\) process, and \(\theta\) is the speed of the mean reversion of \(X\) process towards \(\mu\). Finally, \(dW_t\) is the increment of a standard wiener process at interval \((t, t+dt)\) under a real probability measure \(\mathbb{P}\) and is gaussian in nature with a mean value \(0\) and variance \(t\).

In order to make forecasts, calibration of the model is needed. There are three accepted methods including; the least squares as in (Smith, 2010); the maximum log-likelihood (Van den Berg, 2011); and lastly, the jackknife technique (Phillips, 2005). Based on this, Chaiyapo and Phewchean (2017) calibrated their model with the three methods, whereas (Tanaka & Montero, 2017) used only the least squares method, and Barlow (2002) while modeling electricity prices used the log-likelihood one.

However, despite the consistent parameters, a common limitation with the O-U process is the fact that the spot prices can have negative values. Therefore, inspired by this aspect, Ross (1995), Bessembinder et al. (1995) and Schwartz (1997) proposed a modification of the geometric Ornstein–Uhlenbeck developed previously by Dixit and Pyndick (1994) i.e. the exponential O-U process.

Discussion of O-U Process

Here, we derive the analytical solution to the stochastic differential equation for the Ornstein-Uhlenbeck process: \(\displaystyle dX_{t}=\kappa (\theta -X_{t})\,dt+\sigma \,dW_{t}\) where the \(W_t\) is a standard wiener process, and \(\kappa > 0\), \(\theta \ and \ \sigma > 0\) Motivated by the observation that \(\theta\) is supposed to be the long-term mean of the process \(X_t\), we can simplify the above SDE by introducing the change of variable; \(Y_t = X_t - \theta\) that subtracts off the mean. The \(Y_t\) satisfies the SDE; \(dY_t = dX_t = -\kappa Y_t dt + \sigma dW_t\) for this SDE, the process \(Y_t\) is seen to have a drift towards zero value, at an exponential rate \(\kappa\). This motivates the change of variables; \(Y_t = e^{-\kappa t} Z_t \ \Leftrightarrow \ Z_t = e^{\kappa t} Y_t\) , which should remove the drift. With the application of the product rule for the It\(\tilde{o}\) integral; $dZ_t = e^{t} Y_t dt + e^{t} dY_t $ $= e^{t} Y_t dt + e^{t}(-Y_t dt + dW_t) $ \(= 0 dt + \sigma e^{\kappa t} dW_t\) The solution for \(Z_t\) is then obtained by It\(\tilde{o}\)-intergrating both sides from s to t; \(Z_t = Z_s + \sigma \int_{s}^{t} e^{\kappa u} dW_u\) Reversing the changes of variables, we have; \(Y_t = e^{-\kappa t } Z_t = e^{-\kappa(t-s)} Y_s + \sigma e^{-\kappa t} \int_{s}^{t} e^{\kappa u} dW_u\) and; \(X_t = Y_t + \theta = \theta + e^{-\kappa(t-s)}(X_s - \theta) + \sigma \int_{s}^{t} e^{-\kappa (t-u)} dW_u\)

Therefore, assuming $ X_{s}$ is constant, the mean becomes; \({E} (X_{t})=X_{s}\exp^{-\kappa (t-s)}+\theta (1-\exp^{-\kappa (t-s)})\) Next, using the It\(\tilde{o}\) isometry the covariance function is given by;

\(Cov(X_{s},X_{t})= \frac{\sigma ^{2} }{2\theta } \left(e^{-\theta |t-s|}-e^{-\theta (t+s)}\right)\)

For long-term moments, the mean will tend towards \(\theta\) i.e. the long-term mean, while the variance will be given by \(\frac{\sigma ^{2} }{2\theta }\)

Methodology - Calibration of OU-Process

The choice of the stochastic process to represent the uncertainties, in this case the O-U process, is supported by theoretical considerations referenced in theory and by previous works as discussed in the review with evidence of strong mean-reversion patterns. As stated earlier, the Ornstein-Uhlenbeck (Vasicek process) is stationary, Gaussian, and Markovian. Here, the calibration will be done through the maximum log-likelihood method which works by searching the estimates that maximize the joint loglikelihood function using the first-order conditions.

carbon_em <- read.csv(file.choose(), header = T, sep = ",")
head(carbon_em)

##   date co2_em_kt
## 1 1960  2427.554
## 2 1961  2401.885
## 3 1962  2625.572
## 4 1963  2856.593
## 5 1964  2827.257
## 6 1965  2467.891

Below is a representation of the carbon emissions in Kenya;

plot(carbon_em, type = "l", lwd = 2.5, xlab = "Time", ylab = "Carbon Emission", col = 1, main = "Historical Kenyan Carbon Emission Levels")

Calibration using Maximum Likelihood estimates by employing a function;

ouFit.ML = function(spread) {
    n = length(spread)
    delta = 1  # delta
    Sx = sum(spread[1:n - 1])
    Sy = sum(spread[2:n])
    Sxx = sum((spread[1:n - 1])^2)
    Syy = sum((spread[2:n])^2)
    Sxy = sum((spread[1:n - 1]) * (spread[2:n]))
    mu = (Sy * Sxx - Sx * Sxy)/((n - 1) * (Sxx - Sxy) - (Sx^2 - Sx * Sy))
    theta = -log((Sxy - mu * Sx - mu * Sy + (n - 1) * mu^2)/(Sxx - 2 * mu * 
        Sx + (n - 1) * mu^2))/delta
    a = exp(-theta * delta)
    sigmah2 = (Syy - 2 * a * Sxy + a^2 * Sxx - 2 * mu * (1 - a) * (Sy - a * 
        Sx) + (n - 1) * mu^2 * (1 - a)^2)/(n - 1)
    sigma = sqrt((sigmah2) * 2 * theta/(1 - a^2))
    theta = list(theta = theta, mu = mu, sigma = sigma, sigmah2 = sigmah2)
    return(theta)
}

data <- as.vector(carbon_em[,2])
head(data)

## [1] 2427.554 2401.885 2625.572 2856.593 2827.257 2467.891

ouFit.ML(data)

## $theta
## [1] -0.01220075
## 
## $mu
## [1] -11601.64
## 
## $sigma
## [1] 782.1778
## 
## $sigmah2
## [1] 619327.6

paraData <- c(ouFit.ML(data)[[1]], ouFit.ML(data)[[2]], ouFit.ML(data)[[3]])
paraData

## [1] -1.220075e-02 -1.160164e+04  7.821778e+02

# Simple mean-reverting process:% dx = nu (xbar - x) dt + sigma dz
# generate a function for simulation
OU.sim <- function(para = c(0.01, 200, 1), x0= 1, periods=100){
  #periods=100; #Number of periods
theta = para[1]; # speed of reversion
mu = para[2];
sigma = para[3]; # sigma in monthly terms
# theta = 1
sigma2 = sigma * sqrt((1-exp(-2*theta))/(2*theta));   # dt=1;
x= as.vector(rep(0, periods))  # all zero, 100 vector
epsilon=rnorm(periods)
x[1]= mu; # Starting value of first row of x
for (i in 2:periods){
  x[i] = x[i-1] + mu*(1-exp(-theta)) + x[i-1]*(exp(-theta)-1) + sigma2*epsilon[i-1];
}

return(x) 
}

simdata <- OU.sim(paraData, x0 = 2427.554)

plot(simdata, type = "l", lwd = 2.5, xlab = "Time", ylab = "Carbon Emissions", col = 1, main = "Simulation of Carbon Emission")

We can employ the yuima package also in the simulations. Recall that the package uses this notation: \[ dS_t = \theta(\mu - X_t)dt + \sigma dW_t \]. The delta of t \[ dt=(1/n)=(1/100) \] by default. We can change it by using \(grid=setSampling(Terminal=1,n=1000)\). We would like to generate 1000 random samples with 1000 time periods for carbon emissions using the parameters from the ml-estimates.

require(yuima)

grid = setSampling(Terminal = 1, n = 1000)

## Warning in yuima.warn("'delta' (re)defined."): 
## YUIMA: 'delta' (re)defined.

m1 = setModel(drift = "theta*(mu-x)", diffusion = "sigma", state.var = "x", 
    time.var = "t", solve.var = "x", xinit = ouFit.ML(data)[[2]])
Xdata = simulate(m1, true.param = list(mu = ouFit.ML(data)[[2]], sigma = ouFit.ML(data)[[3]], theta = ouFit.ML(data)[[1]]), sampling = grid)
plot(Xdata, main = "Simulation of Carbon Emission Levels in Kenya")

Findings

From the calibration, we find the unknown parameters by MLE as; theta : -0.01220075 (this represents the speed at which the process is pulled towards the long-term mean), mu : -11601.64 (this represents the long-term mean that the process reverts to), sigma and sigmah2 as 782.1778 and 619327.6 respectively (both represent the extent of volatility of the process, i.e. deviations from the mean). The same parameter estimates are then employed in the simulations above.

Conclusion

This borrowed strategy shows that we can use ouFit.ML function to estimate the parameters for the O-U process \(\theta\), \(\mu\), and \(\sigma\). And them use the “yuima” package to generate the simulative price vector/matrices. The O-U process does a pretty good work on our dataset. Nonetheless, future works might include restrictions to some of the parameters, as some of them cannot have negative values and the present methodology does not impose those conditions. Also, as the maximum log-likelihood method produces a punctual estimator for each parameter, further extensions could be to define intervals for each estimator given a confidence level.

References

https://rstudio-pubs-static.s3.amazonaws.com/19584_ce31e798cffb430982fd2f8979b1a87f.html

https://www.statisticshowto.com/wp-content/uploads/2016/01/Calibrating-the-Ornstein.pdf

M.T., Barlow(2002): A diffusion model for electricity prices. Math. Finance 12(4), 287–298 (2002)

F.E., Benth & K.H., Karlsen(2005): A note on Merton’s portfolio selection problem for the Schwartz mean-reversion model. Stoch. Anal. Appl. 23(4), 687–704

H., Bessembinder etal., (1995): Mean reversion in equilibrium asset prices: evidence from the futures term structure. J. Finance 50(1), 361–375

N., Chaiyapo & N., Phewchean(2017): An application of Ornstein–Uhlenbeck process to commodity pricing in Thailand. Adv.Differ. Equ.

A.K., Dixit & R.S., Pindyck(1994): Investment Under Uncertainty. Princeton University Press, Princeton

D.R., Ribeiro & S.D., Hodges(2004): Equilibrium model for commodity prices: competitive and monopolistic markets. Working Paper

S., Ross(1995): Hedging long run commitments: exercises in incomplete market pricing. Working paper

E.S., Schwartz(1997): The stochastic behavior of commodity prices: implications for valuation and hedging. J. Finance 52(3), 923–973

W., Smith(2010): On the simulation and estimation of mean-reverting Ornstein–Uhlenbeck process: especially as applied to commodities markets and modelling

A.T.,Tanaka & C.M., Carrasco(2017): Valorización de opciones reales: el modelo Ornstein–Uhlenbeck. J. Econ. Finance Adm. Sci. 21(41), 56–62

P.C.B., Phillips & J., Yu(2005): Jackknifing bond option prices. Rev. Financ. Stud. 18(2), 707–742 (2005)

Van den Berg(2011): Calibrating the Ornstein–Uhlenbeck (Vasicek) model.