Этот пример показывает, как моделировать цены на электроэнергию с помощью возвращающейся среднее значение модели с сезонностью и компонентом скачка. Модель калибруется под реальной вероятностью с помощью исторических цен на электроэнергию. Рыночная цена риска получена из цен фьючерсов. Нейтральная к риску симуляция Монте-Карло проводится с помощью калиброванной модели и рыночной цены риска. Результаты симуляции используются, чтобы оценить бермудскую опцию с ценами на электроэнергию как базовое.
Скачки выставки цен на электроэнергию в ценах временами высокого спроса, когда дополнительный, менее эффективных методов производства электроэнергии, принесены онлайн, чтобы обеспечить достаточное электроснабжение. Кроме того, у них есть видный сезонный компонент, наряду с возвращением, чтобы означать уровни. Поэтому эти характеристики должны быть включены в модель цен на электроэнергию [2].
В этом примере цена на электроэнергию моделируется как:
где спотовая цена электричества. Логарифм цены на электроэнергию моделируется с двумя компонентами: и . Компонент детерминированная сезонная часть модели, и стохастическая часть модели. Тригонометрические функции используются к модели можно следующим образом [3]:
где постоянные параметры, и пересчитанные на год факторы времени. Стохастический компонент моделируется как процесс Орнстейна-Ахленбека (среднее возвращение) со скачками:
Параметры и параметры возвращения к среднему уровню. Параметр энергозависимость, и стандартное Броуновское движение. Размер скачка , с нормально распределенным средним значением , и стандартное отклонение . Пуассоновский процесс имеет интенсивность скачка .
Демонстрационные цены на электроэнергию с 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);
Во-первых, детерминированная часть сезонности калибруется с помощью метода наименьших квадратов. Поскольку функция сезонности линейна относительно параметров , оператор наклонной черты влево (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)');
Второй этап должен калибровать стохастическую часть. Модель для потребности, которые будут дискретизированы, чтобы провести калибровку. Чтобы дискретизировать, примите, что существует Бернуллиевый процесс для событий скачка. Таким образом, существует самое большее один скачок в день, поскольку этот пример калибрует против ежедневных цен на электроэнергию. Дискретизированное уравнение:
с вероятностью и,
с вероятностью , где и независимые стандартные нормальные случайные переменные, и . Функция плотности данный [1,4]:
Параметры может быть калиброван путем минимизации отрицательной логарифмической функции правдоподобия:
Первое ограничение неравенства, , эквивалентно . Колебания и mustBePositive. В последнем неравенстве, между 0 и 1, потому что это представляет вероятность скачка, происходящего в время. В этом примере примите это один день. Поэтому за один год существует самое большее 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);
Чтобы калибровать рыночную цену риска против значений фьючерсов рынка, используйте следующее уравнение:
где наблюдаемое значение фьючерсов во время , и ожидаемое значение под реальной мерой во время . Уравнение было получено с помощью той же методологии, как описано в [3]. Этот пример принимает, что рыночная цена риска полностью управляется Броуновским движением. Рыночная цена риска, , может быть решен путем дискретизации вышеупомянутого уравнения и решения системы линейных уравнений.
% 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;
Однажды получен, нейтральная к риску симуляция может быть проведена с помощью следующей динамики:
с вероятностью и
с вероятностью .
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.