В этом примере показано, как симулировать цены на электроэнергию с помощью возвращающейся среднее значение модели с сезонностью и компонентом скачка. Модель калибруется под реальной вероятностью с помощью исторических цен на электроэнергию. Рыночная цена риска получена из цен фьючерсов. Нейтральная к риску симуляция Монте-Карло проводится с помощью калиброванной модели и рыночной цены риска. Результаты симуляции используются, чтобы оценить бермудскую опцию с ценами на электроэнергию как базовое.
Скачки выставки цен на электроэнергию в ценах временами высокого спроса, когда дополнительный, менее эффективных методов производства электроэнергии, принесены онлайн, чтобы обеспечить достаточное электроснабжение. Кроме того, у них есть видный сезонный компонент, наряду с возвращением, чтобы означать уровни. Поэтому эти характеристики должны быть включены в модель цен на электроэнергию [2].
В этом примере цена на электроэнергию моделируется как:
где спотовая цена электричества. Логарифм цены на электроэнергию моделируется с двумя компонентами: и . Компонент детерминированная сезонная часть модели, и стохастическая часть модели. Тригонометрические функции используются к модели можно следующим образом [3]:
где постоянные параметры, и пересчитанные на год факторы времени. Стохастический компонент моделируется как процесс Орнстейна-Ахленбека (среднее возвращение) со скачками:
Параметры и параметры возвращения к среднему уровню. Параметр энергозависимость, и стандартное Броуновское движение. Размер скачка , с нормально распределенным средним значением , и стандартное отклонение . Пуассоновский процесс имеет интенсивность скачка .
Демонстрационные цены на электроэнергию с 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 ($)');
% Obtain the log of prices. logPrices = log(Prices); % Obtain the 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 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)');
Второй этап должен калибровать стохастическую часть. Модель для потребности, которые будут дискретизированы, чтобы провести калибровку. Чтобы дискретизировать, примите, что существует Бернуллиевый процесс для событий скачка. Таким образом, существует самое большее один скачок в день, поскольку этот пример калибрует против ежедневных цен на электроэнергию. Дискретизированное уравнение:
с вероятностью и,
с вероятностью , где и независимые стандартные нормальные случайные переменные, и . Функция плотности данный [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 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.3358
Калиброванные параметры и дискретизированная модель позволяют нам симулировать цены на электроэнергию под реальной вероятностью. Симуляция проводится в течение приблизительно 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');
До этой точки параметры были калиброваны под реальной вероятностью. Однако к ценовым опциям, вам нужна симуляция под нейтральной к риску вероятностью. Чтобы получить это, вычислите рыночную цену риска от цен фьючерсов, чтобы вывести нейтральные к риску параметры. Предположим, что существуют ежемесячные фьючерсные контракты, доступные на рынке, которые улаживаются ежедневно в течение месяца контракта. Например, такие фьючерсы для рынка электроэнергии 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);
Чтобы калибровать рыночную цену риска против значений фьючерсов рынка, используйте следующее уравнение:
где наблюдаемое значение фьючерсов во время , и ожидаемое значение под реальной мерой во время . Уравнение было получено с помощью той же методологии как описано в [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 the 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 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');
Нейтральные к риску симулированные значения используются в качестве входа в функциональный 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.
asianbykv
| asianbylevy
| asianbyls
| asiansensbykv
| asiansensbylevy
| asiansensbyls
| lookbackbycvgsg
| lookbackbyls
| lookbacksensbycvgsg
| lookbacksensbyls
| optpricebysim
| optstockbyblk
| optstockbyls
| optstocksensbyblk
| optstocksensbyls
| spreadbybjs
| spreadbyfd
| spreadbykirk
| spreadbyls
| spreadsensbybjs
| spreadsensbyfd
| spreadsensbykirk
| spreadsensbyls