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

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

Обзор модели

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

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

log(Pt)=f(t)+Xt

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

f(t)=s1sin(2πt)+s2потому что(2πt)+s3sin(4πt)+s4потому что(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 electricity prices and futures prices
load('electricity_prices.mat');

% Plot electricity prices
figure;
plot(PriceDates, Prices);
datetick();
title('Electricity Prices');
xlabel('Date');
ylabel('Price ($)');

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

% Obtain 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 log price and seasonality line
figure;
subplot(2, 1, 1);
plot(PriceDates, logPrices);
datetick();
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);
datetick();
title('log(price) with Seasonality Removed');
xlabel('Date');
ylabel('log(Prices)');

Второй этап должен калибровать стохастическую часть. Модель для 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=1Tжурнал(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 daily prices
dt = 1/365;

% PDF for 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 maximum likelihood
params = mle(Pt,'pdf',mrjpdf,'start',x0,'lowerbound',lb,'upperbound',ub,...
    'optimfun','fmincon');

% Obtain 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 = daysadd(PriceDates(end),0:nPeriods-1);
SimPriceTimes = yearfrac(PriceDates(1), SimPriceDates);
CSim = seasonMatrix(SimPriceTimes);
logSimPrices = SimPrices + repmat(CSim*seasonParam,1,nTrials);

% Plot 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;
datetick();
title('Actual log(price) and Simulated log(price)');
xlabel('Date');
ylabel('log(price)');
legend('market', 'simulation');

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

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

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

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

% Obtain 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 expected futures price under real-world measure
SimPricesExp = mean(PricesSim, 2);

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

журнал(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 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 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');
datetick();
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);
datetick();
title('Market Price of Risk');
xlabel('Date');
ylabel('Market Price of Risk');

Оценка бермудской опции

Нейтральные к риску симулированные значения используются в качестве входа в функциональный optpricebysim в Financial Instruments Toolbox™, чтобы оценить европейца, бермудца или американскую опцию на ценах на электроэнергию. Ниже, цена вычисляется для бермудского колл-опциона 2D года с двумя возможностями осуществления. Первое осуществление после одного года, и второе в зрелости опции.

% Settle, maturity of option
Settle = FutValuationDate;
Maturity = addtodate(FutValuationDate, 2, 'year');

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

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

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

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

Больше о

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

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