exponenta event banner

Моделирование цен на электроэнергию со средней реверсией и скачкообразной диффузией

В этом примере показано, как моделировать цены на электроэнергию с использованием модели возврата среднего значения с сезонностью и компонентом скачка. Модель калибруется под реальной вероятностью с использованием исторических цен на электроэнергию. Рыночная цена риска получается из фьючерсных цен. Нейтральное к риску моделирование Монте-Карло проводится с использованием калиброванной модели и рыночной цены риска. Результаты моделирования используются для определения цены на бермудский вариант с ценами на электроэнергию в качестве базового.

Обзор модели

Цены на электроэнергию резко возрастают в периоды высокого спроса, когда в режиме реального времени вводятся дополнительные, менее эффективные методы выработки электроэнергии для обеспечения достаточного предложения электроэнергии. Кроме того, они имеют заметный сезонный компонент, наряду с обращением к средним уровням. Поэтому эти характеристики должны быть включены в модель цен на электроэнергию [2].

В этом примере цена электроэнергии моделируется следующим образом:

log (Pt) = f (t) + Xt

где Pt - спотовая цена электроэнергии. Логарифм цены электроэнергии моделируется двумя компонентами: f (t) и Xt. Компонент f (t) является детерминированной сезонной частью модели, а Xt - стохастической частью модели. Тригонометрические функции используются для моделирования f (t) следующим образом [3]:

f (t) = s1sin (2.dt) + s2cos (2.dt) + s3sin (4.dt) + s4cos (4.dt) + s5

где si, i = 1,..., 5 - постоянные параметры, а t - ежегодные временные коэффициенты. Стохастический компонент Xt смоделирован как процесс Орнштейна-Уленбека (средний реверс) со скачками:

dXt = (α-αXt) dt + λ dWt + J (мкJ,

Параметры α и λ являются средними параметрами реверсии. Параметр λ - волатильность, а Wt - стандартное броуновское движение. Размер скачка - J (мкДж, в соответствии с J), с нормально распределенным средним мкДж и стандартным отклонением (в соответствии с J). Процесс Пуассона Δ (λ) имеет интенсивность скачка λ.

Цены на электроэнергию

Выборочные цены на электроэнергию с 1 января 2010 года по 11 ноября 2013 года загружены и нанесены на график ниже. Prices содержать цены на электроэнергию, и PriceDates содержат даты, связанные с ценами. Рассчитывается логарифм цен и годовые временные коэффициенты.

% 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. The axes with title Electricity Prices 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);

Калибровка

Во-первых, детерминированная часть сезонности калибруется методом наименьших квадратов. Поскольку функция сезонности является линейной по отношению к параметрам si, оператор обратной косой черты (mldivide) используется. После калибровки сезонность удаляется из логарифма цены. Ниже представлен логарифм тенденций цен и сезонности. Также наносится десезонализированный логарифм цены.

% 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. Axes 1 with title log(price) and Seasonality contains 2 objects of type line. These objects represent log(Price), seasonality. Axes 2 with title log(price) with Seasonality Removed contains an object of type line.

Второй этап - калибровка стохастической части. Модель для Xt должна быть дискретизирована для проведения калибровки. Для дискретизации предположим, что существует процесс Бернулли для событий прыжка. То есть, есть максимум один скачок в день, поскольку этот пример калибруется по ежедневным ценам на электроэнергию. Дискретизированное уравнение:

Xt = αΔt + ϕXt-1 +

с вероятностью (1-λΔt) и,

Xt = αΔt + ϕXt-1 +

с вероятностью λ Δt, где λ и λ J являются независимыми стандартными нормальными случайными переменными, и, ((") = 1-Δt. Функция плотности Xt данного Xt-1 составляет [1,4]:

f (Xt 'Xt-1) = (λ Δt) N1 (Xt' Xt-1) + (1-λΔt) N2 (Xt 'Xt-1)

N1 (Xt'Xt-1) = ((σ2 +σJ2))-12exp (-(Xt \U 03B1\\U 0394\t \U 03D5\Xt 1 \U 03BC\J) 22 (σ2 +σJ2))

N2 (Xt 'Xt-1) = (2íα2) -12exp (- (Xt-αΔt-ϕXt-1) 22start2)

Параметры λ = {α, start, мкJ, start2, σJ2, λ} могут быть откалиброваны путём минимизации функции отрицательного логарифмического правдоподобия:

minθ-∑t=1Tlog (f (Xt 'Xt-1))

subjecttoϕ<1,σ2>0,σJ2>0,0≤λΔt≤1

Первое ограничение неравенства, ϕ <1, эквивалентно κ> 0. Летучие вещества λ и λ J должны быть положительными. В последнем неравенстве λ Δt находится между 0 и 1, поскольку представляет вероятность скачка, происходящего за Δt времени. В этом примере предположим, что Δt равен одному дню. Поэтому есть максимум 365 прыжков за один год. mle функция из Toolbox™ статистики и машинного обучения хорошо подходит для решения вышеуказанной задачи максимального правдоподобия.

% 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

Моделирование Монте-Карло

Калиброванные параметры и дискретизированная модель позволяют нам моделировать цены на электроэнергию под реальной вероятностью. Моделирование проводится в течение приблизительно 2 лет с 10000 испытаний. Он превышает 2 года, чтобы включить все даты в последний месяц моделирования. Это связано с тем, что ожидаемые цены моделирования для даты истечения срока действия фьючерсного контракта необходимы в следующем разделе для расчета рыночной цены риска. Сезонность добавляется обратно к моделируемым путям. График для одной траектории моделирования отображается ниже.

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. Axes 1 with title Actual log(price) and Simulated log(price) contains 3 objects of type line. These objects represent market, simulation. Axes 2 with title Actual Prices and Simulated Prices contains 2 objects of type line. These objects represent market, simulation.

Калибровка рыночной цены риска

До этого момента параметры были откалиброваны под реальной вероятностью. Однако для ценовых вариантов необходимо моделирование под нейтральной для риска вероятностью. Для этого вычислите рыночную цену риска из фьючерсных цен, чтобы получить нейтральные для риска параметры. Предположим, что на рынке имеются ежемесячные фьючерсные контракты, которые рассчитываются ежедневно в течение контрактного месяца. Например, такие фьючерсы на рынок электроэнергии PJM перечислены на Чикагской товарной бирже [5].

Фьючерсы рассчитываются ежедневно в течение месяца контракта. Поэтому можно получить дневную стоимость фьючерсов, предположив, что стоимость фьючерсов является постоянной для месяца контракта. Ожидаемые фьючерсные цены из реального показателя также необходимы для расчета рыночной цены риска. Это можно получить из моделирования, проведенного в предыдущем разделе.

% 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);

Для калибровки рыночной цены риска по рыночной стоимости фьючерсов используйте следующее уравнение:

журнал (FtET) =-σe-kt∫0teksmsds

где Ft - наблюдаемая фьючерсная стоимость в момент времени t, а Et - ожидаемая величина при реальном измерении в момент времени t. Уравнение было получено с использованием той же методологии, что описана в [3]. В этом примере предполагается, что рыночная цена риска полностью зависит от броуновского движения. Рыночная цена риска, mt, может быть решена путем дискретизации вышеуказанного уравнения и решения системы линейных уравнений.

% 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;

Моделирование нейтральных по риску цен

После получения mt можно провести нейтральное с точки зрения риска моделирование с использованием следующей динамики:

Xt = αΔt + ϕXt-1-σmt-1Δt +

с вероятностью (1-λΔt) и

Xt = αΔt + ϕXt-1-σmt-1Δt +

с вероятностью λ Δ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);

Ожидаемые значения из нейтрального к риску моделирования строятся на основе рыночных фьючерсных значений. Это подтверждает, что нейтральное к риску моделирование точно воспроизводит рыночные фьючерсные ценности.

% 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. Axes 1 with title Market Futures Prices and Simulated Futures Prices contains 2 objects of type line. These objects represent market, simulation. Axes 2 with title Market Price of Risk contains an object of type line.

Ценообразование бермудского опциона

В качестве входных данных функции используются нейтральные для риска моделируемые значения. optpricebysim в Toolbox™ «Финансовые инструменты» цена европейского, бермудского или американского варианта цены на электроэнергию. Ниже цена рассчитывается для двухлетнего варианта вызова на Бермудские острова с двумя возможностями упражнений. Первое упражнение проводится через год, а второе - на момент погашения опциона.

% 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

Ссылки

[1] Эскрибано, Альваро, Пена, Хуан Игнасио, Вильяплана, Пабло. «Моделирование цен на электроэнергию: международные доказательства». Университет Карлоеса III, Мадрид, рабочий документ 02-27, 2002 год.

[2] Люсия, Хулио Дж., Шварц, Эдуаро. «Цены на электроэнергию и производные энергии: доказательства от Северной биржи электроэнергии». Обзор деривативных исследований. Том 5, выпуск 1, стр. 5-50, 2002.

[3] Зайферт, Ян, Ухриг-Хомбург, Марлиезе. «Моделирование скачков цен на электроэнергию: теория и эмпирические данные». Обзор деривативных исследований. Том 10, стр. 59-85, 2007.

[4] Вильяплана, Пабло. «Деривативы ценовой силы: двухфакторный подход» скачок-диффузия «». Университет Карлоеса III, Мадрид, рабочий документ 03-18, 2003 год.

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

См. также

| | | | | | | | | | | | | | | | | | | | | |

Связанные примеры

Подробнее

Внешние веб-сайты