- Added Solutions to Chapters 6-11 (together with Python codes)

Solutions to Exercises/Chapter 10/Python Codes/Exercise_10_5.py

@@ -0,0 +1,97 @@
+#%%
+"""
+Created on 29 Nov 2019
+Conditional expectation
+@author: Lech A. Grzelak
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+def GeneratePathsTwoStocksEuler(NoOfPaths,NoOfSteps,T,r,S10,S20,rho,sigma1,sigma2):
+    Z1 = np.random.normal(0.0,1.0,[NoOfPaths,NoOfSteps])
+    Z2 = np.random.normal(0.0,1.0,[NoOfPaths,NoOfSteps])
+    W1 = np.zeros([NoOfPaths, NoOfSteps+1])
+    W2 = np.zeros([NoOfPaths, NoOfSteps+1])
+    # Initialization
+    X1 = np.zeros([NoOfPaths, NoOfSteps+1])
+    X1[:,0] =np.log(S10)
+    X2 = np.zeros([NoOfPaths, NoOfSteps+1])
+    X2[:,0] =np.log(S20)
+    time = np.zeros([NoOfSteps+1])
+    dt = T / float(NoOfSteps)
+
+    for i in range(0,NoOfSteps):
+        # Making sure that samples from a normal have mean 0 and variance 1
+        if NoOfPaths > 1:
+            Z1[:,i] = (Z1[:,i] - np.mean(Z1[:,i])) / np.std(Z1[:,i])
+            Z2[:,i] = (Z2[:,i] - np.mean(Z2[:,i])) / np.std(Z2[:,i])
+            Z2[:,i] = rho *Z1[:,i] + np.sqrt(1.0-rho**2.0)*Z2[:,i]
+        W1[:,i+1] = W1[:,i] + np.power(dt, 0.5)*Z1[:,i]
+        W2[:,i+1] = W2[:,i] + np.power(dt, 0.5)*Z2[:,i]
+        X1[:,i+1] = X1[:,i] + (r -0.5*sigma1**2.0)* dt + sigma1 * (W1[:,i+1] -
+          W1[:,i])
+        X2[:,i+1] = X2[:,i] + (r -0.5*sigma2**2.0)* dt + sigma2 * (W2[:,i+1] -
+              W2[:,i])
+        time[i+1] = time[i] +dt
+
+    # Return stock paths
+    paths = {"time":time,"S1":np.exp(X1),"S2":np.exp(X2)}
+    return paths
+
+def getEVBinMethod(S,v,NoOfBins):
+    if (NoOfBins != 1):
+        mat = np.transpose(np.array([S,v]))
+        # Sort all the rows according to the first column
+        val = mat[mat[:,0].argsort()]
+        binSize = int((np.size(S)/NoOfBins))
+        expectation = np.zeros([np.size(S),2])
+        for i in range(1,binSize-1):
+            sample = val[(i-1)*binSize:i*binSize,1]
+            expectation[(i-1)*binSize:i*binSize,0] = val[(i-1)*binSize:i*binSize,0]
+            expectation[(i-1)*binSize:i*binSize,1] = np.mean(sample)
+    return expectation
+
+def main():
+    NoOfPaths = 10000
+    NoOfSteps = 100
+    T         = 5.0
+    rho       = -0.5
+    y10       = 1.0
+    y20       = 0.05
+
+    # Set 1
+    sigma1    = 0.3
+    sigma2    = 0.3
+
+    # Set 2
+    #sigma1    = 0.9
+    #sigma2    = 0.9
+
+    paths = GeneratePathsTwoStocksEuler(NoOfPaths, NoOfSteps, T, 0.0, y10, y20, rho, sigma1, sigma2)
+    y1 = paths["S1"]
+    y2 = paths["S2"]
+
+    # Analytical expression for the conditional expectation    
+    condE = lambda y1: y20 * (y1/y10)**(rho*sigma2/sigma1)*np.exp(0.5*T*(rho*sigma2*sigma1-sigma2**2*rho**2))
+
+    plt.figure(1)
+    plt.grid()
+    plt.plot(y1[:,-1],y2[:,-1],'.')
+    y1Grid = np.linspace(np.min(y1[:,-1]), np.max(y1[:,-1]), 2500)
+    plt.plot(y1Grid,condE(y1Grid),'r')
+
+    # Bin method
+    E = getEVBinMethod(y1[:,-1], y2[:,-1], 50)
+    plt.plot(E[:,0],E[:,1],'k')
+    plt.legend(['samples','E[Y2|Y1]-Monte Carlo','E[Y2|Y1]-Analytical'])
+    plt.xlabel('Y1')
+    plt.ylabel('Y2')
+
+    # for y1 = 1.75 we have
+    y1 = 1.75
+    print("Analytical expression for y1={0} yields E[Y2|Y1={0}] = {1}".format(y1,condE(y1)))
+
+    condValueInterp = lambda y1: np.interp(y1,E[:,0],E[:,1])
+    print("Monte Carlo, for y1={0} yields E[Y2|Y1={0}] = {1}".format(y1,condValueInterp(y1)))
+
+main()

Solutions to Exercises/Chapter 10/Python Codes/Exercise_10_7.py

@@ -0,0 +1,207 @@
+#%%
+"""
+Created on Jan 8 2020
+Exercise 10.7: The Bates model and pricing of forward start options
+@author: Lech A. Grzelak
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.stats as st
+import enum 
+import scipy.optimize as optimize
+
+# set i= imaginary number
+i   = np.complex(0.0,1.0)
+
+# This class defines puts and calls
+class OptionType(enum.Enum):
+    CALL = 1.0
+    PUT = -1.0
+
+def CallPutOptionPriceCOSMthd_FrwdStart(cf,CP,r,T1,T2,K,N,L):
+    # cf   - characteristic function as a functon, in the book denoted as \varphi
+    # CP   - C for call and P for put
+    # S0   - Initial stock price
+    # r    - interest rate (constant)
+    # K    - list of strikes
+    # N    - Number of expansion terms
+    # L    - size of truncation domain (typ.:L=8 or L=10)  
+
+    tau = T2 - T1
+    # reshape K to a column vector
+    if K is not np.array:
+        K = np.array(K).reshape([len(K),1])
+
+    # Adjust strike
+    K = K + 1.0
+
+    #assigning i=sqrt(-1)
+    i = np.complex(0.0,1.0) 
+    x0 = np.log(1.0 / K)   
+
+    # truncation domain
+    a = 0.0 - L * np.sqrt(tau)
+    b = 0.0 + L * np.sqrt(tau)
+
+    # sumation from k = 0 to k=N-1
+    k = np.linspace(0,N-1,N).reshape([N,1])  
+    u = k * np.pi / (b - a);  
+
+    # Determine coefficients for Put Prices  
+    H_k = CallPutCoefficients(CP,a,b,k)   
+    mat = np.exp(i * np.outer((x0 - a) , u))
+    temp = cf(u) * H_k 
+    temp[0] = 0.5 * temp[0]    
+    value = np.exp(-r * T2) * K * np.real(mat.dot(temp))     
+    return value
+
+# Determine coefficients for Put and Call Prices 
+def CallPutCoefficients(CP,a,b,k):
+    if CP==OptionType.CALL:                  
+        c = 0.0
+        d = b
+        coef = Chi_Psi(a,b,c,d,k)
+        Chi_k = coef["chi"]
+        Psi_k = coef["psi"]
+        if a < b and b < 0.0:
+            H_k = np.zeros([len(k),1])
+        else:
+            H_k      = 2.0 / (b - a) * (Chi_k - Psi_k)  
+    elif CP==OptionType.PUT:
+        c = a
+        d = 0.0
+        coef = Chi_Psi(a,b,c,d,k)
+        Chi_k = coef["chi"]
+        Psi_k = coef["psi"]
+        H_k      = 2.0 / (b - a) * (- Chi_k + Psi_k)               
+
+    return H_k    
+
+def Chi_Psi(a,b,c,d,k):
+    psi = np.sin(k * np.pi * (d - a) / (b - a)) - np.sin(k * np.pi * (c - a)/(b - a))
+    psi[1:] = psi[1:] * (b - a) / (k[1:] * np.pi)
+    psi[0] = d - c
+
+    chi = 1.0 / (1.0 + np.power((k * np.pi / (b - a)) , 2.0)) 
+    expr1 = np.cos(k * np.pi * (d - a)/(b - a)) * np.exp(d)  - np.cos(k * np.pi 
+                  * (c - a) / (b - a)) * np.exp(c)
+    expr2 = k * np.pi / (b - a) * np.sin(k * np.pi * 
+                        (d - a) / (b - a))   - k * np.pi / (b - a) * np.sin(k 
+                        * np.pi * (c - a) / (b - a)) * np.exp(c)
+    chi = chi * (expr1 + expr2)
+
+    value = {"chi":chi,"psi":psi }
+    return value
+
+# Forward start Black-Scholes option price
+def BS_Call_Option_Price_FrwdStart(K,sigma,T1,T2,r):
+    if K is list:
+        K = np.array(K).reshape([len(K),1])
+    K = K + 1.0 
+    tau = T2 - T1
+    d1    = (np.log(1.0 / K) + (r + 0.5 * np.power(sigma,2.0))* tau) / (sigma * np.sqrt(tau))
+    d2    = d1 - sigma * np.sqrt(tau)
+    value = np.exp(-r*T1) * st.norm.cdf(d1) - st.norm.cdf(d2) * K * np.exp(-r * T2)
+    return value
+
+# Implied volatility for the forward start call option
+def ImpliedVolatility_FrwdStart(marketPrice,K,T1,T2,r):
+    # To determine initial volatility we interpolate define a grid for sigma
+    # and interpolate on the inverse
+    sigmaGrid = np.linspace(0,2,200)
+    optPriceGrid = BS_Call_Option_Price_FrwdStart(K,sigmaGrid,T1,T2,r)
+    sigmaInitial = np.interp(marketPrice,optPriceGrid,sigmaGrid)
+    print("Initial volatility = {0}".format(sigmaInitial))
+
+    # Use determined input for the local-search (final tuning)
+    func = lambda sigma: np.power(BS_Call_Option_Price_FrwdStart(K,sigma,T1,T2,r) - marketPrice, 1.0)
+    impliedVol = optimize.newton(func, sigmaInitial, tol=1e-15)
+    print("Final volatility = {0}".format(impliedVol))
+    return impliedVol
+
+def ChFBatesModelForwardStart(r,T1,T2,kappa,gamma,vbar,v0,rho,xiP,muJ,sigmaJ):
+    i = np.complex(0.0,1.0)
+    tau = T2 - T1
+    D1 = lambda u: np.sqrt(np.power(kappa-gamma*rho*i*u,2)+(u*u+i*u)*gamma*gamma)
+    g  = lambda u: (kappa-gamma*rho*i*u-D1(u))/(kappa-gamma*rho*i*u+D1(u))
+    C  = lambda u: (1.0-np.exp(-D1(u)*tau))/(gamma*gamma*(1.0-g(u)*np.exp(-D1(u)*tau)))\
+        *(kappa-gamma*rho*i*u-D1(u))
+    # Note that we exclude the term -r*tau, as the discounting is performed in the COS method
+    AHes  = lambda u: r*i*u*tau + kappa*vbar*tau/gamma/gamma *(kappa-gamma*rho*i*u-D1(u))\
+        - 2*kappa*vbar/gamma/gamma*np.log((1.0-g(u)*np.exp(-D1(u)*tau))/(1.0-g(u)))
+    A = lambda u: AHes(u) - xiP * i * u * tau *(np.exp(muJ+0.5*sigmaJ*sigmaJ) - 1.0) + \
+            xiP * tau * (np.exp(i*u*muJ - 0.5 * sigmaJ * sigmaJ * u * u) - 1.0)
+    c_bar = lambda t1,t2: gamma*gamma/(4.0*kappa) * (1.0 - np.exp(-kappa*(t2-t1)))
+    delta = 4.0*kappa*vbar/gamma/gamma
+    kappa_bar = lambda t1, t2: 4.0*kappa*v0*np.exp(-kappa*(t2-t1))/(gamma*gamma*(1.0-np.exp(-kappa*(t2-t1))))
+    term1 = lambda u: A(u) + C(u)*c_bar(0.0,T1)*kappa_bar(0.0,T1)/(1.0-2.0*C(u)*c_bar(0.0,T1))
+    term2 = lambda u: np.power(1.0/(1.0-2.0*C(u)*c_bar(0.0,T1)),0.5*delta)
+    cf = lambda u: np.exp(term1(u)) * term2(u)
+    return cf
+
+def mainCalculation():
+    CP  = OptionType.CALL
+    r   = 0.00
+
+    TMat1=[[1.0,3.0],[2.0,4.0],[3.0, 5.0],[4.0, 6.0]]
+    TMat2=[[1.0,2.0],[1.0,3.0],[1.0, 4.0],[1.0, 5.0]]
+
+    K = np.linspace(-0.4,4.0,50)
+    K = np.array(K).reshape([len(K),1])
+
+    N = 500
+    L = 10
+
+    # Bates model parameters
+    kappa = 0.6
+    gamma = 0.2
+    vbar  = 0.1
+    rho   = -0.5
+    v0    = 0.05   
+
+    muJ   = 0.05
+    sigmaJ= 0.2
+    xiP   = 0.1
+
+    plt.figure(1)
+    plt.grid()
+    plt.xlabel('strike, K')
+    plt.ylabel('implied volatility')
+    legend = []
+    for T_pair in TMat1:
+       T1= T_pair[0]
+       T2= T_pair[1]
+       cf = ChFBatesModelForwardStart(r,T1,T2,kappa,gamma,vbar,v0,rho,xiP,muJ,sigmaJ)
+       # Forward-start option from the COS method
+       valCOS = CallPutOptionPriceCOSMthd_FrwdStart(cf,CP,r,T1,T2,K,N,L)
+       # Implied volatilities
+       IV =np.zeros([len(K),1])
+       for idx in range(0,len(K)):
+           IV[idx] = ImpliedVolatility_FrwdStart(valCOS[idx],K[idx],T1,T2,r)
+       plt.plot(K,IV*100.0)
+       legend.append('T1={0} & T2={1}'.format(T1,T2))
+    plt.legend(legend)
+
+    plt.figure(2)
+    plt.grid()
+    plt.xlabel('strike, K')
+    plt.ylabel('implied volatility')
+    legend = []
+    for T_pair in TMat2:
+       T1= T_pair[0]
+       T2= T_pair[1]
+       cf = ChFBatesModelForwardStart(r,T1,T2,kappa,gamma,vbar,v0,rho,xiP,muJ,sigmaJ)
+
+       # Forward-start option from the COS method
+       valCOS = CallPutOptionPriceCOSMthd_FrwdStart(cf,CP,r,T1,T2,K,N,L)
+
+       # Implied volatilities
+       IV =np.zeros([len(K),1])
+       for idx in range(0,len(K)):
+           IV[idx] = ImpliedVolatility_FrwdStart(valCOS[idx],K[idx],T1,T2,r)
+       plt.plot(K,IV*100.0)
+       legend.append('T1={0} & T2={1}'.format(T1,T2))
+
+    plt.legend(legend)
+
+mainCalculation()