Simulating Electricity Prices with Mean-Reversion and Jump-Diffusion

Open Live Script

This example shows how to simulate electricity prices using a mean-reverting model with seasonality and a jump component. The model is calibrated under the real-world probability using historical electricity prices. The market price of risk is obtained from futures prices. A risk-neutral Monte Carlo simulation is conducted using the calibrated model and the market price of risk. The simulation results are used to price a Bermudan option with electricity prices as the underlying.

Overview of the Model

Electricity prices exhibit jumps in prices at periods of high demand when additional, less efficient electricity generation methods, are brought on-line to provide a sufficient supply of electricity. In addition, they have a prominent seasonal component, along with reversion to mean levels. Therefore, these characteristics should be incorporated into a model of electricity prices [2].

In this example, electricity price is modeled as:

$l o g (P_{t}) = f (t) + X_{t}$

where $P_{t}$ is the spot price of electricity. The logarithm of electricity price is modeled with two components: $f (t)$ and $X_{t}$ . The component $f (t)$ is the deterministic seasonal part of the model, and $X_{t}$ is the stochastic part of the model. Trigonometric functions are used to model $f (t)$ as follows [3]:

$f (t) = s_{1} \sin (2 π t) + s_{2} \cos (2 π t) + s_{3} \sin (4 π t) + s_{4} \cos (4 π t) + s_{5}$

where $s_{i}, i = 1, . . ., 5$ are constant parameters, and $t$ is the annualized time factors. The stochastic component $X_{t}$ is modeled as an Ornstein-Uhlenbeck process (mean-reverting) with jumps:

$d X_{t} = (α - κ X_{t}) d t + σ d W_{t} + J (μ_{J}, σ_{J}) d Π (λ)$

The parameters $α$ and $κ$ are the mean-reversion parameters. Parameter $σ$ is the volatility, and $W_{t}$ is a standard Brownian motion. The jump size is $J (μ_{J}, σ_{J})$ , with a normally distributed mean $μ_{J}$ , and a standard deviation $σ_{J}$ . The Poisson process $Π (λ)$ has a jump intensity of $λ$ .

Electricity Prices

Sample electricity prices from January 1, 2010 to November 11, 2013 are loaded and plotted below. Prices contain the electricity prices, and PriceDates contain the dates associated with the prices. The logarithm of the prices and annual time factors are calculated.

% Load the electricity prices and futures prices.
load('electricity_prices.mat');
PriceDates = datetime(PriceDates,'ConvertFrom','datenum');
FutExpiry = datetime(FutExpiry,'ConvertFrom','datenum');
FutValuationDate = datetime(FutValuationDate,'ConvertFrom','datenum');
% Plot the electricity prices.
figure;
plot(PriceDates, Prices);
title('Electricity Prices');
xlabel('Date');
ylabel('Price ($)');

Figure contains an axes object. The axes object with title Electricity Prices, xlabel Date, ylabel Price ($) contains an object of type line.

% Obtain the log of prices.
logPrices = log(Prices);

% Obtain the annual time factors from dates.
PriceTimes = yearfrac(PriceDates(1), PriceDates);

Calibration

First, the deterministic seasonality part is calibrated using the least squares method. Since the seasonality function is linear with respect to the parameters $s_{i}$ , the backslash operator (mldivide) is used. After the calibration, the seasonality is removed from the logarithm of price. The logarithm of price and seasonality trends are plotted below. Also, the de-seasonalized logarithm of price is plotted.

% Calibrate parameters for the seasonality model.
seasonMatrix = @(t) [sin(2.*pi.*t) cos(2.*pi.*t) sin(4.*pi.*t) ...
    cos(4.*pi.*t) t ones(size(t, 1), 1)];
C = seasonMatrix(PriceTimes);
seasonParam = C\logPrices;

% Plot the log price and seasonality line.
figure;
subplot(2, 1, 1);
plot(PriceDates, logPrices);
title('log(price) and Seasonality');
xlabel('Date');
ylabel('log(Prices)');
hold on;
plot(PriceDates, C*seasonParam, 'r');
hold off;
legend('log(Price)', 'seasonality');

% Plot de-seasonalized log price
X = logPrices-C*seasonParam;
subplot(2, 1, 2);
plot(PriceDates, X);
title('log(price) with Seasonality Removed');
xlabel('Date');
ylabel('log(Prices)');

Figure contains 2 axes objects. Axes object 1 with title log(price) and Seasonality, xlabel Date, ylabel log(Prices) contains 2 objects of type line. These objects represent log(Price), seasonality. Axes object 2 with title log(price) with Seasonality Removed, xlabel Date, ylabel log(Prices) contains an object of type line.

The second stage is to calibrate the stochastic part. The model for $X_{t}$ needs to be discretized to conduct the calibration. To discretize, assume that there is a Bernoulli process for the jump events. That is, there is at most one jump per day since this example is calibrating against daily electricity prices. The discretized equation is:

$X_{t} = α Δ t + ϕ X_{t - 1} + σ ξ$

with probability $(1 - λ Δ t)$ and,

$X_{t} = α Δ t + ϕ X_{t - 1} + σ ξ + μ_{J} + σ_{J} ξ_{J}$

with probability $λ Δ t$ , where $ξ$ and $ξ_{J}$ are independent standard normal random variables, and $ϕ = 1 - κ Δ t$ . The density function of $X_{t}$ given $X_{t - 1}$ is [1,4]:

$f (X_{t} | X_{t - 1}) = (λ Δ t) N_{1} (X_{t} | X_{t - 1}) + (1 - λ Δ t) N_{2} (X_{t} | X_{t - 1})$

$N_{1} (X_{t} | X_{t - 1}) = (2 π (σ^{2} + σ_{J}^{2}))^{- \frac{1}{2}} \exp (\frac{- (X_{t} - α Δ t - ϕ X_{t - 1} - μ_{J})^{2}}{2 (σ^{2} + σ_{J}^{2})})$

$N_{2} (X_{t} | X_{t - 1}) = (2 π σ^{2})^{- \frac{1}{2}} \exp (\frac{- (X_{t} - α Δ t - ϕ X_{t - 1})^{2}}{2 σ^{2}})$

The parameters $θ = {α, ϕ, μ_{J}, σ^{2}, σ_{J}^{2}, λ}$ can be calibrated by minimizing the negative log likelihood function:

$m i n_{θ} - \sum_{t = 1}^{T} \log (f (X_{t} | X_{t - 1}))$

$s u b j e c t t o ϕ < 1, σ^{2} > 0, σ_{J}^{2} > 0, 0 \leq λ Δ t \leq 1$

The first inequality constraint, $ϕ < 1$ , is equivalent to $κ > 0$ . The volatilities $σ$ and $σ_{J}$ must be positive. In the last inequality, $λ Δ t$ is between 0 and 1, because it represents the probability of a jump occurring in $Δ t$ time. In this example, assume that $Δ t$ is one day. Therefore, there is at most 365 jumps in one year. The mle function from the Statistics and Machine Learning Toolbox™ is well suited to solve the above maximum likelihood problem.

% Prices at t, X(t).
Pt = X(2:end);

% Prices at t-1, X(t-1).
Pt_1 = X(1:end-1);

% Discretization for the daily prices.
dt = 1/365;

% PDF for the discretized model.
mrjpdf = @(Pt, a, phi, mu_J, sigmaSq, sigmaSq_J, lambda) ...
    lambda.*exp((-(Pt-a-phi.*Pt_1-mu_J).^2)./ ...
    (2.*(sigmaSq+sigmaSq_J))).* (1/sqrt(2.*pi.*(sigmaSq+sigmaSq_J))) + ...
    (1-lambda).*exp((-(Pt-a-phi.*Pt_1).^2)/(2.*sigmaSq)).* ...
    (1/sqrt(2.*pi.*sigmaSq));

% Constraints: 
% phi < 1 (k > 0)
% sigmaSq > 0
% sigmaSq_J > 0
% 0 <= lambda <= 1
lb = [-Inf -Inf -Inf 0 0 0];
ub = [Inf 1 Inf Inf Inf 1];

% Initial values.
x0 = [0 0 0 var(X) var(X) 0.5];

% Solve the maximum likelihood.
params = mle(Pt,'pdf',mrjpdf,'start',x0,'lowerbound',lb,'upperbound',ub,...
    'optimfun','fmincon');

% Obtain the calibrated parameters.
alpha = params(1)/dt

alpha = 
-20.1060

kappa = (1-params(2))/dt

kappa = 
188.2535

mu_J = params(3)

mu_J = 
0.2044

sigma = sqrt(params(4)/dt);
sigma_J = sqrt(params(5))

sigma_J = 
0.2659

lambda = params(6)/dt

lambda = 
98.3357

Monte Carlo Simulation

The calibrated parameters and the discretized model allow us to simulate electricity prices under the real-world probability. The simulation is conducted for approximately 2 years with 10,000 trials. It exceeds 2 years to include all the dates in the last month of simulation. This is because the expected simulation prices for the futures contract expiry date is required in the next section to calculate the market price of risk. The seasonality is added back on the simulated paths. A plot for a single simulation path is plotted below.

rng default;

% Simulate for about 2 years.
nPeriods = 365*2+20;
nTrials = 10000;
n1 = randn(nPeriods,nTrials);
n2 = randn(nPeriods, nTrials);
j = binornd(1, lambda*dt, nPeriods, nTrials);
SimPrices = zeros(nPeriods, nTrials);
SimPrices(1,:) = X(end);
for i=2:nPeriods
    SimPrices(i,:) = alpha*dt + (1-kappa*dt)*SimPrices(i-1,:) + ...
                sigma*sqrt(dt)*n1(i,:) + j(i,:).*(mu_J+sigma_J*n2(i,:));
end

% Add back seasonality.
SimPriceDates = PriceDates(end) + days(0:(nPeriods-1))';
SimPriceTimes = yearfrac(PriceDates(1), SimPriceDates);
CSim = seasonMatrix(SimPriceTimes);
logSimPrices = SimPrices + repmat(CSim*seasonParam,1,nTrials);

% Plot the logarithm of Prices and simulated logarithm of Prices.
figure;
subplot(2, 1, 1);
plot(PriceDates, logPrices);
hold on;
plot(SimPriceDates(2:end), logSimPrices(2:end,1), 'red');
seasonLine = seasonMatrix([PriceTimes; SimPriceTimes(2:end)])*seasonParam;
plot([PriceDates; SimPriceDates(2:end)], seasonLine, 'green');
hold off;
title('Actual log(price) and Simulated log(price)');
xlabel('Date');
ylabel('log(price)');
legend('market', 'simulation');

% Plot the prices and simulated prices.
PricesSim = exp(logSimPrices);
subplot(2, 1, 2);
plot(PriceDates, Prices);
hold on;
plot(SimPriceDates, PricesSim(:,1), 'red');
hold off;
title('Actual Prices and Simulated Prices');
xlabel('Date');
ylabel('Price ($)');
legend('market', 'simulation');

Figure contains 2 axes objects. Axes object 1 with title Actual log(price) and Simulated log(price), xlabel Date, ylabel log(price) contains 3 objects of type line. These objects represent market, simulation. Axes object 2 with title Actual Prices and Simulated Prices, xlabel Date, ylabel Price ($) contains 2 objects of type line. These objects represent market, simulation.

Calibration of the Market Price of Risk

Up to this point, the parameters were calibrated under the real-world probability. However, to price options, you need the simulation under the risk-neutral probability. To obtain this, calculate the market price of risk from futures prices to derive the risk-neutral parameters. Suppose that there are monthly futures contracts available on the market, which are settled daily during the contract month. For example, such futures for the PJM electricity market are listed on the Chicago Mercantile Exchange [5].

The futures are settled daily during the contract month. Therefore, you can obtain daily futures values by assuming the futures value is constant for the contract month. The expected futures prices from the real-world measure are also needed to calculate the market price of risk. This can be obtained from the simulation conducted in the previous section.

% Obtain the daily futures prices.
FutPricesDaily = zeros(size(SimPriceDates));
for i=1:nPeriods
    idx = find(year(SimPriceDates(i)) == year(FutExpiry) & ...
        month(SimPriceDates(i)) == month(FutExpiry));
    FutPricesDaily(i) = FutPrices(idx);
end

% Calculate the expected futures price under real-world measure.
SimPricesExp = mean(PricesSim, 2);

To calibrate the market price of risk against market futures values, use the following equation:

$\log (\frac{F_{t}}{E_{t}}) = - σ e^{- k t} \int_{0}^{t} e^{k s} m_{s} d s$

where $F_{t}$ is the observed futures value at time $t$ , and $E_{t}$ is the expected value under the real-world measure at time $t$ . The equation was obtained using the same methodology as described in [3]. This example assumes that the market price of risk is fully driven by the Brownian motion. The market price of risk, $m_{t}$ , can be solved by discretizing the above equation and solving a system of linear equations.

% Setup system of equations.
t0 = yearfrac(PriceDates(1), FutValuationDate);
tz = SimPriceTimes-t0;
b = -log(FutPricesDaily(2:end) ./ SimPricesExp(2:end)) ./ ...
    (sigma.*exp(-kappa.*tz(2:end)));
A = (1/kappa).*(exp(kappa.*tz(2:end)) - exp(kappa.*tz(1:end-1)));
A = tril(repmat(A', size(A,1), 1));

% Precondition to stabilize numerical inversion.
P = diag(1./diag(A));
b = P*b;
A = P*A;

% Solve for the market price of risk.
riskPremium = A\b;

Simulation of Risk-Neutral Prices

Once $m_{t}$ is obtained, risk-neutral simulation can be conducted using the following dynamics:

$X_{t} = α Δ t + ϕ X_{t - 1} - σ m_{t - 1} Δ t + σ ξ$

with probability $(1 - λ Δ t)$ and

$X_{t} = α Δ t + ϕ X_{t - 1} - σ m_{t - 1} Δ t + σ ξ + μ_{J} + σ_{J} ξ_{J}$

with probability $λ Δ t$ .

nTrials = 10000;
n1 = randn(nPeriods, nTrials);
n2 = randn(nPeriods, nTrials);
j = binornd(1, lambda*dt, nPeriods, nTrials);

SimPrices = zeros(nPeriods, nTrials);
SimPrices(1,:) = X(end);
for i=2:nPeriods
    SimPrices(i,:) = alpha*dt + (1-kappa*dt)*SimPrices(i-1,:) + ...
        sigma*sqrt(dt)*n1(i,:) - sigma*dt*riskPremium(i-1) + ...
        j(i,:).*(mu_J+sigma_J*n2(i,:));
end

% Add back seasonality.
CSim = seasonMatrix(SimPriceTimes);
logSimPrices = SimPrices + repmat(CSim*seasonParam,1,nTrials);

% Convert the log(Price) to Price.
PricesSim = exp(logSimPrices);

The expected values from the risk-neutral simulation are plotted against the market futures values. This confirms that the risk-neutral simulation closely reproduces the market futures values.

% Obtain expected values from the risk-neutral simulation.
SimPricesExp = mean(PricesSim,2);
fexp = zeros(size(FutExpiry));
for i = 1:size(FutExpiry,1)
    idx = SimPriceDates == FutExpiry(i);    
    if sum(idx)==1
        fexp(i) = SimPricesExp(idx);
    end
end

% Plot expected values from the simulation against market futures prices.
figure;
subplot(2,1,1);
plot(FutExpiry, FutPrices(1:size(FutExpiry,1)),'-*');
hold on;
plot(FutExpiry, fexp, '*r');
hold off;
title('Market Futures Prices and Simulated Futures Prices');
xlabel('Date');
ylabel('Price');
legend('market', 'simulation', 'Location', 'NorthWest');
subplot(2,1,2);
plot(SimPriceDates(2:end), riskPremium);
title('Market Price of Risk');
xlabel('Date');
ylabel('Market Price of Risk');

Figure contains 2 axes objects. Axes object 1 with title Market Futures Prices and Simulated Futures Prices, xlabel Date, ylabel Price contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent market, simulation. Axes object 2 with title Market Price of Risk, xlabel Date, ylabel Market Price of Risk contains an object of type line.

Pricing a Bermudan Option

The risk-neutral simulated values are used as input into the function optpricebysim in the Financial Instruments Toolbox™ to price a European, Bermudan, or American option on electricity prices. Below, the price is calculated for a two-year Bermudan call option with two exercise opportunities. The first exercise is after one year, and the second is at the maturity of the option.

% Settle, maturity of option.
Settle = FutValuationDate;
Maturity = FutValuationDate + calyears(2);

% Create the interest-rate term structure.
riskFreeRate = 0.01;
Basis = 0;
Compounding = -1;
RateSpec = intenvset('ValuationDate', Settle, 'StartDates', Settle, ...
    'EndDates', Maturity, 'Rate', riskFreeRate, 'Compounding', ...
    Compounding, 'Basis', Basis);

% Cutoff the simulation at maturity.
endIdx = find(SimPriceDates == Maturity);
SimPrices = PricesSim(1:endIdx,:);
Times = SimPriceTimes(1:endIdx) - SimPriceTimes(1);

% Bermudan call option with strike 60, two exercise opportunities, after
% one year and at maturity.
OptSpec = 'call';
Strike = 60;
ExerciseTimes = [Times(366) Times(end)];
Price = optpricebysim(RateSpec, SimPrices, Times, OptSpec, Strike, ...
    ExerciseTimes)

Price = 
1.1085

References

[1] Escribano, Alvaro, Pena, Juan Ignacio, Villaplana, Pablo. "Modeling Electricity Prices: International Evidence." Universidad Carloes III de Madrid, Working Paper 02-27, 2002.

[2] Lucia, Julio J., Schwartz, Eduaro. "Electricity Prices and Power Derivatives: Evidence from the Nordic Power Exchange." Review of Derivatives Research. Vol. 5, Issue 1, pp 5-50, 2002.

[3] Seifert, Jan, Uhrig-Homburg, Marliese. "Modelling Jumps in Electricity Prices: Theory and Empirical Evidence." Review of Derivatives Research. Vol. 10, pp 59-85, 2007.

[4] Villaplana, Pablo. "Pricing Power Derivatives: A Two-Factor Jump-Diffusion Approach." Universidad Carloes III de Madrid, Working Paper 03-18, 2003.

[5] https://www.cmegroup.com