# 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.

• Thank you @noob2, the problem is how to pass from those terminal conditions to the boundary conditions to apply when executing the numerical method? In my lecture notes it is written that the right (hence assuming $S>K$) terminal condition for a European call option at each step $n$ is $V_t = S_t-K\exp(-r\,m\,\Delta t)$ (where $m$ is the time step) where the term $\exp(-r\,m\,\Delta t)$ is, if I have understood well, the discounting factor (factor by which a future cash flow must be multiplied in order to obtain the present value)...continue in next comment – sound wave Sep 18 '20 at 15:41
• @noob2 In my notes it is also written that the forward contract with price at delivery $F(S,T)$, where $S$ is the current price at time 0 and $T$ is the time to expiry, is equivalent to having an option with $\text{payoff}(S_T) = S_T-K$ with $K=F(S,T)$ (I think $K$ is called like strike price). From what I understood what you wrote $V_T=S_T$ is the right terminal condition for the forward contract, which is different from the condition of the option equivalent to the forward contract. Applying the discounting factor we obtain the same formula as in the European call, but with a different $K$? – sound wave Sep 18 '20 at 15:47
• I am very sorry, my mistake, for a forward $\text{payoff}(S_T)=S_T−K$ is correct. – noob2 Sep 18 '20 at 16:13
• For a forward, you can have an analytical formula. Why do you still need the finite difference approach? – Gordon Sep 21 '20 at 12:46
• @Gordon I was asked to solve this in order to test the goodness of the finite difference method (and then to compare it also with a monte carlo simulation) – sound wave Sep 21 '20 at 13:46