Matlab codes for Solving financial differential equations using differentiation matrices

Americandemo
This script solves the American put option pricing problem as in section 3.

bsdemo
This script solves the European call option pricing problem as in section 2.

blackscholes
This function computes the Black-Scholes formula for a call option. The code returns an array for array-valued S and t arguments. The values at expiry (t=T) are computed using the payoff formula.

get_ode_blackscholes
This function computes the terminal condition and coefficient matrices and vectors for the european call option DAE in section 2.

b_blackscholes
This function computes the forcing term b(t) for the european call option DAE.

get_ode_heston
This function computes the terminal condition and coefficient matrices and vectors for the stochastic volatility option DAE in section 4.

b_heston
This computes the forcing term b(t) for the stochastic volatility option DAE.

get_diffmatrix
This function computes the first and second order differentiation matrices for finite differences and for spectral collocation. In the case of finite differences, central difference formulas are used for inner nodes and one-sided difference formulas are used for boundary nodes. In the case of spectral collocation, the function is just a front end to chebdif from the Differentiation Matrix Suite, with renumbering to give the nodes in increasing order.

evalsc
This function evaluates the polynomial that interpolates data at the Chebyshev nodes. The function is just a front end to chebint from the Differentiation Matrix Suite, with renumbering to give the nodes in increasing order.

ivpsolver_fi
This function is the fully implicit (backward Euler) solver for a DAE of the form My'=Ay+b(t). The linear system at each time step is solved using direct LU matrix factorization; the sparse factorization is used when the coefficient matrices are sparse, and the factorization is recomputed only when the time step changes.

ivpsolver_cn
Like ivpsolver_fi, but using the Crank-Nicolson (trapezoid) method.

ivpsolver_alex
Like ivpsolver_fi, but using Alexander's implicit Runge-Kutta method.

danglwirl
This function solves the optimal maintenance and shutdown problem from section 5.

Figure_WorkAccuracy_DW
This script produces a Work-Accuracy plot (Figure 5) comparing different discretisations in the solution of the optimal maintenance and shutdown problem.

Figure_WorkAccuracy_BS
This function produces a Work-Accuracy plot (Figure 2) comparing different timestepping methods in the solution of the european call option problem.

Figure_WorkAccuracy_H
This script produces a Work-Accuracy plot (Figure 4) comparing different discretisations in the solution of the stochastic volatility option pricing problem.

heston
This function computes Heston's formula for the value of a call option. The integral is computed with quadl with an error tolerance of 10-8. The singularity of the integrand at the origin is avoided by using the midpoint formula. The infinite upper bound of the integration interval is replaced by 100.

Figure_stability_BS
This script plots errors of numerical solutions of the european call pricing problem (Figure 1).