Explicit Euler Pricing for a Mean-Reverting Commodity Forward

# Solve the Schwartz mean reverting PDE for option pricing using Euler explicit method (matlab)


# Solve the Schwartz mean reverting PDE for option pricing using Euler explicit method (matlab)












> Objective: Implement the Euler Explict Method for solving the PDE for option prices under the Schwartz mean reverting model.

The price evolution of a commodity can be described by the Schwartz SDE $$dS = \alpha(\mu-\log S)Sdt + \sigma S dW, \qquad \begin{array}WW = \text{ Standard Brownian motion} \\ \alpha = \text{ strength of mean reversion}\end{array}$$

Consider the forward contract with price at delivery $F(S,T)$ where $T$ is the expiry time and $S$ is the spot price at time 0 of the underlying commodity following the Schwartz SDE. Such a contract is equivalent to having an option with $\text{payoff}(S_T)=S_T-K$ where $K=F(S,T)$ is the strike price. The PDE for the value of this option can be derived from the Schwartz SDE and its formula is $$\frac{\partial V}{\partial t} + \alpha\Big(\mu-\frac{\mu-r}\alpha -\log S\Big)S\frac{\partial V}{\partial S}+\frac12\sigma^2S^2\frac{\partial^2V}{\partial S^2}-rV = 0$$

A solution of the option value PDE is $V(S,t)=e^{-r(T-t)}F(S,t)$.

Using Euler explicit method, i.e. forward difference on $\dfrac{\partial V}{\partial t}$ and central difference on $\dfrac{\partial V}{\partial S}$ and $\dfrac{\partial^2V}{\partial S^2}$, we obtain an equation of the form (complete formula below) $$\tag{1}V^{n+1}_i = a_i V^n_{i-1} + b_i V^n_i + c_i V^n_{i+1}$$

A simpler case is the Black-Scholes PDE for a European call option, in this case it's easy to compute $S_i$ since they increase linearly ($S_i = i \Delta S$) and the boundary conditions for the Euler method are $V_1 = 0$ and $V_n = S_n-K\exp(-r\ m\ dt)$ (where $m$ is the time step).

But what about the Schwartz PDE? In this case, the option is no more a European call but a forward contract and I don't know how to properly fix all the details: how to compute $S_i$? (are $S_i$ obtained by making the difference between a value of the spot price and the previous one?) Which are the boundary conditions?

This is the MatLab code, it works in the sense that runs without errors but the $V$ approximated with the explicit Euler method is bad since has some huge values ($10^{90}$)

```
% real data of monthly oil spot prices
spot_prices = [ 22.93 15.45 12.61 12.84 15.38 13.43 11.58 15.10 14.87 14.90 15.22 16.11 18.65 17.75 18.30 18.68 19.44 20.07 21.34 20.31 19.53 19.86 18.85 17.27 17.13 16.80 16.20 17.86 17.42 16.53 15.50 15.52 14.54 13.77 14.14 16.38 18.02 17.94 19.48 21.07 20.12 20.05 19.78 18.58 19.59 20.10 19.86 21.10 22.86 22.11 20.39 18.43 18.20 16.70 18.45 27.31 33.51 36.04 32.33 27.28 25.23 20.48 19.90 20.83 21.23 20.19 21.40 21.69 21.89 23.23 22.46 19.50 18.79 19.01 18.92 20.23 20.98 22.38 21.78 21.34 21.88 21.69 20.34 19.41 19.03 20.09 20.32 20.25 19.95 19.09 17.89 18.01 17.50 18.15 16.61 14.51 15.03 14.78 14.68 16.42 17.89 19.06 19.65 18.38 17.45 17.72 18.07 17.16 18.04 18.57 18.54 19.90 19.74 18.45 17.33 18.02 18.23 17.43 17.99 19.03 18.85 19.09 21.33 23.50 21.17 20.42 21.30 21.90 23.97 24.88 23.71 25.23 25.13 22.18 20.97 19.70 20.82 19.26 19.66 19.95 19.80 21.33 20.19 18.33 16.72 16.06 15.12 15.35 14.91 13.72 14.17 13.47 15.03 14.46 13.00 11.35 12.51 12.01 14.68 17.31 17.72 17.92 20.10 21.28 23.80 22.69 25.00 26.10 27.26 29.37 29.84 25.72 28.79 31.82 29.70 31.26 33.88 33.11 34.42 28.44 29.59 29.61 27.24 27.49 28.63 27.60 26.42 27.37 26.20 22.17 19.64 19.39 19.71 20.72 24.53 26.18 27.04 25.52 26.97 28.39 ];

r = .1;    % yearly instantaneous interest rate
T = 1/2;   % expiry time
ts = 3000; % number of time steps
dt = T/ts; % delta t
dS = diff(spot_prices(1:end-1)); % delta S, contains also negative values, is it ok?

alpha = 0.0692;
sigma = 0.0876; % values estimated from data
mu = 3.0582;

S = spot_prices;
S1 = S(2:end-1);     % all but endpoints
t = linspace(0,T,numel(S));
tau = T-t;
lambdahat = (mu-r)/alpha;
muhat = mu-sigma^2/2/alpha-lambdahat;
F = exp( exp(-alpha*tau).*log(S) + muhat*(1-exp(-alpha*tau)) + sigma^2/4/alpha*(1-exp(-2*alpha*tau)) ); % value of the Forward contract
K = F(end); % strike price
V = F .* exp( -r*tau ); value of the option equivalent to the forward contract
plot(t,S)
hold on
plot(t,V)
% Euler explicit method
for m = 1:ts
    V(2:end-1) = .5*dt./dS.*S1 .* ( sigma^2*S1./dS - alpha*mu + (mu-r) + alpha*log(S1) ) .* V(1:end-2) +...
                                                    ( 1 - r*dt - sigma^2*dt./dS.*S1.^2 ) .* V(2:end-1) +...
                 .5*dt./dS.*S1 .* ( sigma^2*S1./dS + alpha*mu - (mu-r) - alpha*log(S1) ) .* V(3:end);
    % boundary conditions
    V(1) = 0; % since there is no money exchanged when signing this contract
    V(end) = S(end)-K*exp(-r*m*dt);
end
plot(t,V)
legend('Spot prices','Option prices','Euler approximation')
```

The full formula for (1) which appears in the for cycle of the code

\begin{align*} V^{n+1}_i &= \frac12\frac{\Delta t}{\Delta S}S_i \Big( \sigma^2\frac{S_i}{\Delta S} -\alpha\mu + (\mu-r) + \alpha \log S_i \Big) V^n_{i-1} \\ &+ \Big(1-r\Delta t - \sigma^2 \frac{\Delta t}{\Delta S^2} S^2_i \Big) V^n_i \\ &+ \frac12\frac{\Delta t}{\Delta S}S_i \Big( \sigma^2\frac{S_i}{\Delta S} +\alpha\mu - (\mu-r) - \alpha \log S_i \Big) V^n_{i+1} \end{align*}

The formula for the price of the forward contract appearing in the code

$$F(S,\tau)=\exp\Big(e^{-\alpha\tau}\log S +(\mu-\frac{\sigma^2}{2\alpha}-\frac{\mu-r}{\alpha})(1-e^{-\alpha\tau})+\frac{\sigma^2}{4\alpha}(1-e^{-2\alpha\tau})\Big)$$

where $\tau=T-t$ is the time to expiry.