Number of Paths for $J$ in Lecture 2 Code

[aanchitnayak]

Hi Prof. Lech,

In the following code:

def GeneratePaths(NoOfPaths,NoOfSteps,S0,T,muJ,sigmaJ,r):    
    # Create empty matrices for Poisson process and for compensated Poisson process
    X = np.zeros([NoOfPaths, NoOfSteps+1])
    S = np.zeros([NoOfPaths, NoOfSteps+1])
    time = np.zeros([NoOfSteps+1])
                
    dt = T / float(NoOfSteps)
    X[:,0] = np.log(S0)
    S[:,0] = S0
    
    Z = np.random.normal(0.0,1.0,[NoOfPaths,NoOfSteps])
    J = np.random.normal(muJ,sigmaJ,[NoOfPaths,NoOfSteps])
    for i in range(0,NoOfSteps):
        # making sure that samples from normal have mean 0 and variance 1
        if NoOfPaths > 1:
            Z[:,i] = (Z[:,i] - np.mean(Z[:,i])) / np.std(Z[:,i])
            
        X[:,i+1]  = X[:,i] + (r - 0.5*J[:,i]**2.0)*dt +J[:,i]*np.sqrt(dt)* Z[:,i]
        time[i+1] = time[i] +dt
        
    S = np.exp(X)
    paths = {"time":time,"X":X,"S":S,"J":J}
    return paths

$J$ is being drawn randomly at each point in the path, making the simulated stochastic process path dependent. In the lecture and the conditional expectation code, we can either select one $j \in J \sim \mathcal{N}(\mu_J, \sigma^2_j)$ for the entire path, or we will have to change the code for the conditional expectation's calculation to add an effective volatility parameter which would integrate the volatility over the path of the price process.

We can either:

Option A: Change the GeneratePaths function

From

    J = np.random.normal(muJ,sigmaJ,[NoOfPaths,NoOfSteps])
    for i in range(0,NoOfSteps):
        # making sure that samples from normal have mean 0 and variance 1
        if NoOfPaths > 1:
            Z[:,i] = (Z[:,i] - np.mean(Z[:,i])) / np.std(Z[:,i])
            
        X[:,i+1]  = X[:,i] + (r - 0.5*J[:,i]**2.0)*dt +J[:,i]*np.sqrt(dt)* Z[:,i]
        time[i+1] = time[i] +dt

to

    J = np.random.normal(muJ,sigmaJ,[NoOfPaths,1])
    for i in range(0,NoOfSteps):
        # making sure that samples from normal have mean 0 and variance 1
        if NoOfPaths > 1:
            Z[:,i] = (Z[:,i] - np.mean(Z[:,i])) / np.std(Z[:,i])
        J_i = J[:, 0]    
        X[:,i+1]  = X[:,i] + (r - 0.5*J_i**2.0)*dt +J_i*np.sqrt(dt)* Z[:,i]
        time[i+1] = time[i] +dt

This would correspond to the idea that while we have stochastic volatility, it is one instance of that volatility over the entire realization of the price process.

Option B: Change the Conditional Expectation Estimate

From

def CallOption_CondExpectation(NoOfPaths,T,S0,K,J,r):
    
    # Jumps at time T
    J_i = J[:,-1]
    
    result = np.zeros([NoOfPaths])
    
    for j in range(0,NoOfPaths):
        sigma = J_i[j]
        result[j] = BS_Call_Put_Option_Price(OptionType.CALL,S0,[K],sigma,0.0,T,r)
        
    return np.mean(result)

to

def call_option_conditional_expectation_timevarying_J(CP, S0, K, T, r, J, dt):
    # J shape: (no_of_paths, no_of_steps)

    V = np.sum(J**2, axis=1) * dt          # integrated variance per path
    sigma_eff = np.sqrt(V / T)             # effective BS vol per path (>= 0)

    # vectorized BS (same as scalar BS but without Python loop)
    K_ = float(np.array(K).reshape(-1)[0])  # ensure scalar strike

    # handle near-zero sigma to avoid divide-by-zero
    eps = 1e-12
    sigma = np.maximum(sigma_eff, eps)

    d1 = (np.log(S0 / K_) + (r + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
    d2 = d1 - sigma*np.sqrt(T)

    if CP == OptionType.CALL:
        prices = st.norm.cdf(d1)*S0 - st.norm.cdf(d2)*K_*np.exp(-r*T)
        # for sigma_eff extremely small, use deterministic limit exactly
        near_zero = sigma_eff < eps
        if np.any(near_zero):
            fwd = S0*np.exp(r*T)
            prices[near_zero] = np.exp(-r*T) * np.maximum(fwd - K_, 0.0)
    elif CP == OptionType.PUT:
        prices = st.norm.cdf(-d2)*K_*np.exp(-r*T) - st.norm.cdf(-d1)*S0
        # for sigma_eff extremely small, use deterministic limit exactly
        near_zero = sigma_eff < eps
        if np.any(near_zero):
            fwd = S0*np.exp(r*T)
            prices[near_zero] = np.exp(-r*T) * np.maximum(K_ - fwd, 0.0)
    else:
        raise ValueError("Unknown CP value")

    return np.mean(prices)

This function would make the estimation of conditional expectation an integration over the path of the stock price process. I am unsure if this is absolutely correct but the code - as it is now - doesn't seem to represent a consistent depiction of the price process as you describe in the video lecture.