## Root finding routines:
##
source('http://waxworksmath.com/Authors/A_F/Conte/Code/Chapter3/utils.R')
##source('../../../../Conte/Code/Chapter3/utils.R')

Nprime = function(x){
    return(exp(-x^2/2)/sqrt(2*pi))
}

N_scalar = function(x){
    ##
    ## x = scalar input (works on only scalar inputs)
    ## 
    gamma = 0.2316419
    k = 1/(1 + gamma*x)
    
    a1 = 0.319381530
    a2 = -0.35653782
    a3 = 1.781477937
    a4 = -1.821255978
    a5 = 1.330274429
    
    if( x < 0 ){
        return(1 - N(-x))
    }else{
        return(1 - Nprime(x)*(a1*k + a2*k^2 + a3*k^3 + a4*k^4 + a5*k^5))
    }
}

N = function(x){
    ##
    ## x = vector inputs
    ##
    return(sapply(x, N_scalar))
}


BS_Call = function(S, t, r, K, T, sigma){
    d1 = ( log(S/K) + (r + sigma^2/2) * (T-t) ) / ( sigma * sqrt(T-t) )
    d2 = ( log(S/K) + (r - sigma^2/2) * (T-t) ) / ( sigma * sqrt(T-t) )
    C = S * N(d1) - K * exp(-r*(T-t)) * N(d2)
    return(C)
}

if( FALSE ){
    ## Verify results in Example 13.2
    ##
    C = BS_Call(62, 0, 0.1, 60, 5/12, 0.2)
    print(round(C, 3))
}


BS_Call_Implied_Volatility = function(C, S0, t, r, K, T){
    sigma_small = 0.01
    sigma_larger = 1.00
    res = bisection_method(sigma_small, sigma_larger, function(x) (BS_Call(S0, t, r, K, T, x) - C), xtol=1.e-6, ftol=le-3)
    return(res)
}


trinomial_lattice_probabilities = function(u, mu, sigma, deltaT){
    ##
    ## This solves the system of linear equations given by Equations 13.23.
    ##
    A = matrix(data=c(1, 1, 1, u, 1, 1/u, u^2, 1, 1/u^2), nrow=3, ncol=3, byrow=T)
    b = c(1, 1 + mu*deltaT, sigma^2*deltaT + (1 + mu*deltaT)^2)
    ps = solve(A, b)
    return(ps)
}