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

Открыть Live Script

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

Обзор модели

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

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

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

где $P_{t}$ - спотовая цена электроэнергии. Логарифм цены на электроэнергию моделируется двумя компонентами: $f (t)$ и $X_{t}$ . Компонент $f (t)$ является детерминированной сезонной частью модели, и $X_{t}$ является стохастической частью модели. Тригонометрические функции используются для моделирования $f (t)$ следующим образом [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}$

где $s_{i}, i = 1, . . ., 5$ являются постоянными параметрами, и $t$ - ежегодные временные факторы. Стохастический компонент $X_{t}$ моделируется как процесс Орнштейна-Уленбека (среднее значение) с переходами:

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

Параметры $α$ и $κ$ являются средними параметрами реверсии. Параметр $σ$ является волатильностью, и $W_{t}$ является стандартным броуновским движением. Размер перехода $J (μ_{J}, σ_{J})$ , с нормально распределенным средним $μ_{J}$ , и стандартное отклонение $σ_{J}$ . The Пуассоновского процесса $Π (λ)$ имеет интенсивность перехода $λ$ .

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

Выборочные цены на электроэнергию с 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);

Калибровка

Во-первых, детерминированная часть сезонности калибруется с помощью метода наименьших квадратов. Поскольку функция сезонности линейна относительно параметров $s_{i}$ , оператор обратной косой черты (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.

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

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

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

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

с вероятностью $λ Δ t$ , где $ξ$ и $ξ_{J}$ являются независимыми стандартными нормальными случайными переменными, и $ϕ = 1 - κ Δ t$ . Функция плотности $X_{t}$ данный $X_{t - 1}$ [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}})$

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

$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$

Первое ограничение неравенства, $ϕ < 1$ , эквивалентно $κ > 0$ . Колебания $σ$ и $σ_{J}$ должен быть положительным. В последнем неравенстве, $λ Δ t$ находится между 0 и 1, потому что это представляет вероятность перехода, происходящего в $Δ t$ время. В этом примере предположим, что $Δ t$ это один день. Поэтому за один год наблюдается самое большее 365 переходы. The 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 (\frac{F_{t}}{E_{t}}) = - σ e^{- k t} \int_{0}^{t} e^{k s} m_{s} d s$

где $F_{t}$ - наблюдаемое значение фьючерса в момент времени $t$ , и $E_{t}$ - ожидаемое значение под действительной мерой в то время $t$ . Уравнение получали по той же методологии, что и в [3]. Этот пример принимает, что рыночная цена риска полностью определяется броуновским движением. Рыночная цена риска, $m_{t}$ , может быть решено путем дискретизации вышеуказанного уравнения и решения системы линейных уравнений.

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

Симуляция нейтральных к риску цен

Однажды $m_{t}$ Получают, нейтральную к риску симуляция можно проводить с помощью следующей динамики:

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

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

$X_{t} = α Δ t + ϕ X_{t - 1} - σ m_{t - 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™ оценить европейскую бермудскую или американскую опцию по ценам на электроэнергию. Ниже, цена рассчитывается для двухлетнего Бермудского вызова опции с двумя возможностями упражнений. Первое упражнение проходит через один год, а второе - по зрелости опции.

% 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] Эскрибано, Альваро, Пена, Хуан Игнасио, Виллаплана, Пабло. Моделирование цен на электроэнергию: международные доказательства. Universidad Carloes III de Madrid, рабочий документ 02-27, 2002.

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

[3] Seifert, Jan, Uhrig-Homburg, Marliese. Моделирование переходов цен на электроэнергию: теория и эмпирические данные. Обзор исследований производных. Том 10, стр 59-85, 2007.

[4] Виллаплана, Пабло. «Ценообразование производных степени: Двухфакторный подход скачка-диффузии». Universidad Carloes III de Madrid, рабочий документ 03-18, 2003.

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

См. также

Подробнее о

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

Торговля энергией и управление рисками с MATLAB (47 мин 31 сек)

Документация