www.gusucode.com > fininst 案例源码程序 matlab代码 > fininst/SimulateElectricityPricesExample.m
%% Simulating Electricity Prices with Mean-Reversion and Jump-Diffusion % % This example shows how to simulate electricity prices using a % mean-reverting model with seasonality and a jump component. The model is % calibrated under the real-world probability using historical % electricity prices. The market price of risk is obtained from futures % prices. A risk-neutral Monte Carlo simulation is conducted using the % calibrated model and the market price of risk. The simulation results are % used to price a Bermudan option with electricity prices as the % underlying. % % Copyright 2014 The MathWorks, Inc. %% Overview of the Model % % Electricity prices exhibit jumps in prices at periods of high demand when % additional, less efficient electricity generation methods are brought % on-line to provide a sufficient supply of electricity. In addition, % they have a prominent seasonal component, along with reversion to mean % levels. Therefore, these characteristics should be incorporated into a % model of electricity prices [2]. % % In this example, electricity price is modeled as: % % $log(P_t) = f(t) + X_t$ % % where $P_t$ is the spot price of electricity. The logarithm of electricity % price is modeled with two components: $f(t)$ and $X_t$. The component % $f(t)$ is the deterministic seasonal part of the model, and $X_t$ is % the stochastic part of the model. Trigonometric functions are used to % model $f(t)$ as follows [3]: % % $f(t) = s_1 \sin(2 \pi t) + s_2 \cos(2 \pi t) + s_3 \sin(4 \pi t) + s_4 % \cos(4 \pi t) + s_5$ % % where $s_i, i=1,...,5$ are constant parameters, and $t$ is the annualized % time factors. The stochastic component $X_t$ is modeled as an % Ornstein-Uhlenbeck process (mean-reverting) with jumps: % % $dX_t = (\alpha - \kappa X_t)dt + \sigma dW_t + J(\mu_J, \sigma_J) % d\Pi(\lambda)$ % % The parameters $\alpha$ and $\kappa$ are the mean-reversion parameters. % Parameter $\sigma$ is the volatility, and $W_t$ is a standard Brownian % motion. The jump size is $J(\mu_J, \sigma_J)$, with normally distributed % mean $\mu_J$ and standard deviation $\sigma_J$. The Poisson process % $\Pi(\lambda)$ has a jump intensity of $\lambda$. %% Electricity Prices % % Sample electricity prices from January 1, 2010 to November 11, 2013 are % loaded and plotted below. |Prices| contain the electricity prices, and % |PriceDates| contain the dates associated with the prices. The logarithm % of the prices and annual time factors are calculated. % 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); %% Calibration % First, the deterministic seasonality part is calibrated using the least % squares method. Since the seasonality function is linear with respect to % the parameters $s_i$, the backslash operator (|mldivide|) is used. After % the calibration, the seasonality is removed from the logarithm of price. % The logarithm of price and seasonality trends are plotted below. Also, % the de-seasonalized logarithm of price is plotted. % 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)'); %% % The second stage is to calibrate the stochastic part. The model for $X_t$ % needs to be discretized in order to conduct the calibration. To % discretize, we assume a Bernoulli process for the jump events. That is, % there is at most one jump per day since we are calibrating against daily % electricity prices. The discretized equation is: % % $X_t=\alpha \Delta t + \phi X_{t-1} + \sigma \xi$ % % with probability $(1 - \lambda \Delta t)$ and, % % $X_t=\alpha \Delta t + \phi X_{t-1} + \sigma \xi + \mu_J + \sigma_J \xi_J$ % % with probability $\lambda \Delta t$, where $\xi$ and $\xi_J$ are % independent standard normal random variables, and $\phi = 1 - \kappa % \Delta t$. The density function of $X_t$ given $X_{t-1}$ is [1,4]: % % $f(X_t|X_{t-1}) = (\lambda \Delta t) N_1(X_t|X_{t-1}) + % (1 - \lambda \Delta t) N_2(X_t|X_{t-1})$ % % $N_1(X_t|X_{t-1}) = (2 \pi (\sigma^2 + \sigma_J^2))^{-\frac{1}{2}} % \exp(\frac{-(X_t - \alpha \Delta t - \phi X_{t-1} - \mu_J)^2}{2 (\sigma^2 % + \sigma_J^2)})$ % % $N_2(X_t|X_{t-1}) = (2 \pi \sigma^2)^{-\frac{1}{2}} \exp(\frac{-(X_t - % \alpha \Delta t - \phi X_{t-1})^2}{2 \sigma^2})$ % % The parameters $\theta = % \{\alpha,\phi,\mu_J,\sigma^2, \sigma_J^2,\lambda\}$ can be calibrated by % minimizing the negative log likelihood function: % % $min_{\theta} % -\sum^T_{t=1} \log(f(X_t|X_{t-1}))$ % % $subject ~to ~\phi < 1, \sigma^2 > 0, \sigma_J^2 > 0, 0 \leq \lambda % \Delta t \leq 1$ % % The first inequality constraint, $\phi < 1$, is equivalent to % $\kappa > 0$. The volatilities $\sigma$ and $\sigma_J$ must be positive. % In the last inequality, $\lambda \Delta t$ is between 0 and 1, because it % represents the probability of a jump occurring in $\Delta t$ time. In this % example, we take $\Delta t$ to be one day. Consequently, there is at most % 365 jumps in one year. The function |mle| from the Statistics and Machine % Learning Toolbox(TM) is well suited to solve the above maximum likelihood % problem. % 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 kappa = params(2)/dt mu_J = params(3) sigma = sqrt(params(4)/dt) sigma_J = sqrt(params(5)) lambda = params(6)/dt %% Monte Carlo Simulation % % The calibrated parameters and the discretized model allow us to simulate % electricity prices under the real-world probability. The simulation is % conducted for approximately 2 years with 10,000 trials. It exceeds 2 % years to include all the dates in the last month of simulation. This is % because the expected simulation prices for the futures contract expiry % date is required in the next section to calculate the market price of % risk. The seasonality is added back on the simulated paths. A plot for % a single simulation path is plotted below. 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'); %% Calibration of the Market Price of Risk % % Up to this point, the parameters were calibrated under the real-world % probability. However, in order to price options, we need the simulation % under the risk-neutral probability. To obtain this, we calculate the % market price of risk from futures prices to derive the risk-neutral % parameters. Suppose that there are monthly futures contracts available on % the market, which are settled daily during the contract month. For % example, such futures for the PJM electricity market is listed on the % Chicago Mercantile Exchange [5]. % % The futures are settled daily during the contract month. Therefore, we % obtain daily futures values by assuming the futures value is constant for % the contract month. The expected futures prices from the real-world % measure are also needed to calculate the market price of risk. This can % be obtained from the simulation conducted in the previous section. % 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); %% % To calibrate the market price of risk against market futures values, we % use the following equation: % % $\log(\frac{F_t}{E_t}) = -\sigma e^{-kt} \int_0^t e^{ks} m_s ds$ % % where $F_t$ is the observed futures value at time $t$, and $E_t$ is the % expected value under the real-world measure at time $t$. The equation was % obtained using the same methodology as described in [3]. We assume that % the market price of risk is fully driven by the Brownian motion. The % market price of risk, $m_t$, can be solved by discretizing the above % equation and solving a system of linear equations. % 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; %% Simulation of Risk-neutral Prices % Once $m_t$ is obtained, risk-neutral simulation can be conducted using % the following dynamics: % % $X_t = \alpha \Delta t + \phi X_{t-1} - \sigma m_{t-1} \Delta t + % \sigma \xi$ % % with probability $(1 - \lambda \Delta t)$ and % % $X_t = \alpha \Delta t + \phi X_{t-1} - \sigma m_{t-1} \Delta t + \sigma % \xi + \mu_J + \sigma_J \xi_J$ % % with probability $\lambda \Delta 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); %% % The expected values from the risk-neutral simulation are plotted against % the market futures values. This confirms that the risk-neutral simulation % closely reproduces the market futures values. % 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'); %% Pricing a Bermudan Option % The risk-neutral simulated values can be used as input into the function % |optpricebysim| in the Financial Instruments Toolbox(TM) to price an % European, Bermudan, or American option on electricity prices. Below, we % price a two year Bermudan call option with two exercise opportunities. % The first excercise is after one year, and the second is at the maturity % of the option. % 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) %% References % % [1] Escribano, Alvaro, Pena, Juan Ignacio, Villaplana, Pablo, Modeling % Electricity Prices: International Evidence, Universidad Carloes III % de Madrid, Working Paper 02-27, 2002. % % [2] Lucia, Julio J., Schwartz, Eduaro, Electricity Prices and Power % Derivatives: Evidence from the Nordic Power Exchange, Review of % Derivatives Research, Vol. 5, Issue 1, pp 5-50, 2002. % % [3] Seifert, Jan, Uhrig-Homburg, Marliese, Modelling Jumps in Electricity % Prices: Theory and Empirical Evidence, Review of Derivatives % Research, Vol. 10, pp 59-85, 2007. % % [4] Villaplana, Pablo, Pricing Power Derivatives: A Two-Factor % Jump-Diffusion Approach, Universidad Carloes III de Madrid, Working % Paper 03-18, 2003. % % [5] http://www.cmegroup.com %