In pedroguarderas/CFINI: Computational finance routines

Introduction

1. Loading packages

library( CFINI )
library( plotly )

A general representation of the equation related to the Black-Scholes models. \begin{equation} \left{ \begin{array}{ll} \dfrac{\partial V}{\partial t}( t, S ) + \alpha( t, S ) \dfrac{\partial^2 V}{\partial S^2}( t, S ) + \beta( t, S ) \dfrac{\partial V}{\partial S}( t, S ) + \gamma( t, S ) V( t, S ) = 0 & \forall ( t, S ) \in [0,T] \times [ S_l, S_h ] \ \text{FC:} & \ V( T, S ) = V_0( S ) & \forall S \in [ S_l, S_h ] \ \text{BC:} & \ V( t, S_l ) = v_1( t ) & \forall t \in [ 0, T ] \ V( t, S_h ) = v_2( t ) & \forall t \in [ 0, T ] \end{array} \right. \end{equation}

the particular case of Black-Scholes is given by the following parameters $\alpha( t, S ) = \frac{1}{2}\sigma^2 S^2$, $\beta( t, S ) = r S$ and $\gamma( t, S ) = -r$ where $\sigma, r$ are constants.

It is important to observe that this model considers a final condition FC, that usually represents the value of the option at the end of the contract.

The idea of being able to produce a differential model that will be capable to produce a consistent pricing model that will hedge the final option is supported by the the second fundamental theorem of asset pricing, which precisely relates the ability to hedge arbitrary claims, to the uniqueness of martingale measures.

Transformation of the particular case of Black-Scholes

The previous model can be transformed in the classical diffusion problem, related to the heat equation, by the right change of variables. \begin{equation} V( t, S ) = e^{\alpha x + \beta \tau} u( x, t )\ \alpha = -\frac{1}{2}\left( \frac{2r}{\sigma^2} - 1 \right) \ \beta = -\frac{1}{4}\left( \frac{2r}{\sigma^2} + 1 \right)^2 \ S = e^x \ t = T - \frac{2 \tau}{\sigma^2} \end{equation}

1. Related coefficients

# Diffusion parameter constant
Nt <- 300
Nx <- 150
alpha <- matrix( 10^(-2.3), Nt, Nx )
theta <- 0.25

1. Time grid

t0 <- 0
t1 <- 1
t <- cf_uniform_grid( t0, t1, Nt )

1. Value grid

x0 <- -0.5
x1 <- 0.5
x <- cf_uniform_grid( x0, x1, Nx )

1. Initial conditions

If <- function( x ) if ( abs(x) <= 0.1 ) return( 1 ) else return( 0 )
If <- Vectorize( If )
I <- sapply( x, FUN = If )

1. Boundary conditions

A <- rep( 0, Nt )
B <- rep( 0, Nt )

\begin{equation} A = \begin{bmatrix} b_1 & a_1 & 0 & 0 & \cdots & 0 \ c_1 & b_2 & a_2 & 0 & \cdots & 0 \ 0 & c_2 & b_3 & a_3 & \cdots & 0 \ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \ 0 & 0 & \cdots & c_{n-2} & b_{n-1} & a_{n-1} \ 0 & 0 & \cdots & 0 & c_{n-1} & b_n \end{bmatrix} \end{equation}

1. Euler implicit scheme has the following form for the equation \begin{equation} u_{n+1,i} - u_{n,i} = \lambda_{n, i} ( u_{n+1,i+1} - 2 u_{n+1,i} + u_{n+1,i-1} ) + \rho_{n,i} ( u_{n+1,i+1/2} - u_{n+1,i-1/2} ) + \gamma_{n,i} u_{n+1,i} \end{equation}

From previous scheme at every time step $n$ we formulate a tridiagonal problem $A u_n = d$, with the folowing definitions. \begin{eqnarray} \lambda_{n, i} & = & \alpha_{n, i} \frac{\Delta t_n}{\Delta x_{i}\ \Delta x_{i+1}} \ \rho_{n, i} & = & \beta_{n, i} \frac{\Delta t_n}{2(x_{i+1} - x_{i})} \ a_i & = & -\lambda_{n,i} \ b_i & = & 1 + 2 \lambda_{n,i} \ c_i & = & -\lambda_{n,i} \ d_i & = & u_{n,i} \end{eqnarray}

Ueu <- cf_diff_solv_euls( alpha, I, A, B, t, x, FALSE )

1. Solving with Crank-Nicolson implicit method \begin{equation} u_{n+1,i} + \theta ( \lambda^2_{n,i} \Delta_x u_{n+1,i+1} - \lambda^1_{n,i} \Delta_x u_{n+1,i} ) = u_{n,i} + ( 1 - \theta ) ( \lambda^2_{n,i} \Delta_x u_{n,i+1} - \lambda^1_{n,i} \Delta_x u_{n,i} ) \end{equation}

\begin{eqnarray} \lambda^1_{n,i} & = & \alpha_{n, i} \frac{\Delta t_n}{\Delta x_{i}\ \Delta x_{i+1}} \ \lambda^2_{n,i} & = & \alpha_{n, i} \frac{\Delta t_n}{\Delta x_{i+1}\ \Delta x_{i+1}} \ a_i & = & -\theta \lambda^1_{n,i} \ b_i & = & 1 + \theta ( \lambda^1_{n,i} + \lambda^2_{n,i} ) \ c_i & = & -\theta \lambda^2_{n,i} \ d_i & = & u_{n,i} + ( 1 - \theta ) ( \lambda^2_{n,i} \Delta_x u_{n,i+1} - \lambda^1_{n,i} \Delta_x u_{n,i} ) \end{eqnarray}

Ucn <- cf_diff_solv_cns( theta, alpha, I, A, B, t, x, FALSE )

alph <- max( alpha )
phi <- function( x ) dnorm( x, mean = 0, sd = sqrt( 2 * alph * t[ Nt ] ) )
s <- sapply( x, FUN = function( y ) integrate( function( x, y ) If( y - x ) * phi( x ), -Inf, Inf, y )$value )

eeu <- matrix( Ueu$u[Nt,] - s, Nx, 1 ) # Euler error
ecn <- matrix( Ucn$u[Nt,] - s, Nx, 1 )

1. Plotting results

plot_ly( x = Ucn$t, y = Ucn$x, z = Ueu$u, alpha = 0.8 ) %>% 
  layout( scene = list( xaxis = list(title = "t"),
                        yaxis = list(title = "S"),
                        zaxis = list(title = "u") ) ) %>%
  add_surface()

pedroguarderas/CFINI documentation built on Feb. 16, 2024, 2:17 p.m.