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

Скрипт Open 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}$ . Пуассоновский процесс $Π (λ)$ имеет интенсивность скачка $λ$ .

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

Демонстрационные цены на электроэнергию с 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) используется. После калибровки сезонность удалена из логарифма цены. Логарифм цены и трендов сезонности построен ниже. Кроме того, 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.

Второй этап должен калибровать стохастическую часть. Модель для $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}$ 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 (\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™, чтобы оценить европейца, бермудца или американскую опцию на ценах на электроэнергию. Ниже, цена вычисляется для бермудского колл-опциона 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

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

Больше о

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

Энергетическая Торговля & управление рисками с MATLAB (47 секунд min 31)

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