Симуляция цен на электроэнергию с возвращением к среднему уровню и диффузией скачка

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

Обзор модели

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

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

log(Pt)=f(t)+Xt

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

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

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

dXt=(α-κXt)dt+σdWt+J(μJ,σJ)dΠ(λ)

Параметры α и κ параметры возвращения к среднему уровню. Параметр σ энергозависимость, и Wt стандартное Броуновское движение. Размер скачка J(μJ,σ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) используется. После калибровки сезонность удалена из логарифма цены. Логарифм цены и трендов сезонности построен ниже. Кроме того, de-seasonalized логарифм цены построен.

% 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+σξ+μJ+σJξJ

с вероятностью λΔ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π(σ2+σJ2))-12exp(-(Xt-αΔt-ϕXt-1-μJ)22(σ2+σJ2))

N2(Xt|Xt-1)=(2πσ2)-12exp(-(Xt-αΔt-ϕXt-1)22σ2)

Параметры θ={α,ϕ,μJ,σ2,σJ2,λ} может быть калиброван путем минимизации отрицательной логарифмической функции правдоподобия:

minθ-t=1Tlog(f(Xt|Xt-1))

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

Первое ограничение неравенства, ϕ<1, эквивалентно κ>0. Колебания σ и σJ mustBePositive. В последнем неравенстве, λΔt между 0 и 1, потому что это представляет вероятность скачка, происходящего в Δt время. В этом примере примите это Δt один день. Поэтому за один год существует самое большее 365 скачков. mle функция от Statistics and Machine Learning 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 лет с 10 000 испытаний. Это превышает 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);

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

log(FtEt)=-σe-kt0teksmsds

где 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+σξ+μJ+σJξJ

с вероятностью λΔ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 в Financial Instruments Toolbox™, чтобы оценить европейца, бермудца или американскую опцию на ценах на электроэнергию. Ниже, цена вычисляется для бермудского колл-опциона 2D года с двумя возможностями осуществления. Первое осуществление после одного года, и второе в зрелости опции.

% 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] Escribano, Альваро, Pena, Хуан Игнасио, Villaplana, Пабло. "Моделирование Цен на электроэнергию: Международное Доказательство". Юниверсидад Карло III де Мадрид, Рабочий документ 02-27, 2002.

[2] Люсия, Хулио Х., Шварц, Eduaro. "Цены на электроэнергию и Производные Степени: Доказательство от скандинавского Exchange Степени". Анализ Исследования Производных. Издание 5, Выпуск 1, стр 5-50, 2002.

[3] Зайферт, январь, Uhrig-хомбург, Marliese. "Моделируя Скачки в Ценах на электроэнергию: Теория и Эмпирическое доказательство". Анализ Исследования Производных. Издание 10, стр 59-85, 2007.

[4] Villaplana, Пабло. "Производные Ценообразования: 2D Факторный Подход Диффузии Скачка". Юниверсидад Карло III де Мадрид, Рабочий документ 03-18, 2003.

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

Смотрите также

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

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

Больше о

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

Для просмотра документации необходимо авторизоваться на сайте