<jupyter_start><jupyter_text>Define function for the computation of the bounds using the penalization approach<jupyter_code>def compute_bound_3dim(payoff,
                       generate_samples_mu1,
                       generate_samples_mu2,
                       generate_samples_mu3,
                        dynamic_options = True,
                        lower_bound=True,
                        batch_size=2**8,
                        gamma = 100,
                        max_magnitude = 1000,
                        max_iter = 1000,
                        nr_neurons = 128,
                        l_r = 0.0001,
                        option_strike1 = 98,
                        option_strike2 = 98,
                        depth = 3,
                        dynamic_factor = 1,
                        lower_epsilon = 0,
                        upper_epsilon = 0):
    # Create Tensors for the Input
    def build_model(dynamic_options):
        x_1 = keras.Input(shape=(1,),name = "x_1")
        x_2 = keras.Input(shape=(1,),name = "x_2")
        x_3 = keras.Input(shape=(1,),name = "x_3")
        x_12 = layers.concatenate([x_1, x_2])
        
        # Create the NN       
        u_1 = layers.Dense(nr_neurons,activation = "relu")(x_1)
        u_2 = layers.Dense(nr_neurons,activation = "relu")(x_2)
        u_3 = layers.Dense(nr_neurons,activation = "relu")(x_3)
        H_1 = layers.Dense(nr_neurons,activation = "relu")(x_1)
        H_2 = layers.Dense(nr_neurons,activation = "relu")(x_12)
        if dynamic_options:
            H_121 = layers.Dense(nr_neurons,activation = "relu")(x_1)
            H_131 = layers.Dense(nr_neurons,activation = "relu")(x_1)
            H_231 = layers.Dense(nr_neurons,activation = "relu")(x_12)
        
        # Create deep layers
        for i in range(depth):
            u_1 = layers.Dense(nr_neurons,activation = "relu")(u_1)
            u_2 = layers.Dense(nr_neurons,activation = "relu")(u_2)
            u_3 = layers.Dense(nr_neurons,activation = "relu")(u_3)
            H_1 = layers.Dense(nr_neurons,activation = "relu")(H_1)            
            H_2 = layers.Dense(nr_neurons,activation = "relu")(H_2)
            if dynamic_options:
                H_121 = layers.Dense(nr_neurons,activation = "relu")(H_121)
                H_131 = layers.Dense(nr_neurons,activation = "relu")(H_131)
                H_231 = layers.Dense(nr_neurons,activation = "relu")(H_231)
          
        # Output Layers
        u_1_out = layers.Dense(1,name = "u_1")(u_1)
        u_2_out = layers.Dense(1,name = "u_2")(u_2)
        u_3_out = layers.Dense(1,name = "u_3")(u_3)
        H_1_out = layers.Dense(1,name = "H_1")(H_1)
        H_2_out = layers.Dense(1,name = "H_2")(H_2)
        if dynamic_options:
            H_121_out = layers.Dense(1,name = "H_121")(H_121)
            H_131_out = layers.Dense(1,name = "H_131")(H_131)
            H_231_out = layers.Dense(1,name = "H_231")(H_231)        
        # Build the model
            model = keras.Model(inputs=[x_1,x_2,x_3],
                                outputs = [u_1_out,u_2_out,u_3_out,H_1_out,H_2_out,
                                           H_121_out,H_131_out,H_231_out])
        else:
            model = keras.Model(inputs=[x_1,x_2,x_3],
                    outputs = [u_1_out,u_2_out,u_3_out,H_1_out,H_2_out])
        return model
    
    # Loss function
    def loss(model,x_1,x_2,x_3,epoch):
        u_1, u_2, u_3, H_1, H_2 = model({"x_1":x_1,"x_2":x_2,"x_3":x_3})
        hedge = u_1 + u_2 + u_3 + H_1*(x_2-x_1) + H_2*(x_3-x_2)
        s = payoff(x_1,x_2,x_3)*(1-2*int(lower_bound)) - hedge
        loss = tf.reduce_mean(u_1+u_2+u_3)+ \
        0.5*gamma_dynamic(gamma, epoch)*tf.reduce_mean(tf.pow(tf.nn.relu(s),2))
        return loss
    
    def loss_dynamic(model,x_1,x_2,x_3,p_121,p_131,p_231,epoch):
        u_1, u_2, u_3, H_1, H_2, H_121, H_131, H_231 = model({"x_1":x_1,"x_2":x_2,"x_3":x_3})
        hedge = u_1 + u_2 + u_3 + H_1*(x_2-x_1) + H_2*(x_3-x_2)
        hedge = hedge + H_121*(tf.nn.relu(x_2-option_strike1)-p_121)+ \
                        H_131*(tf.nn.relu(x_3-option_strike2)-p_131)+\
                        H_231*(tf.nn.relu(x_3-option_strike2)-p_231)
        s = payoff(x_1,x_2,x_3)*(1-2*int(lower_bound)) - hedge
        loss = tf.reduce_mean(u_1+u_2+u_3)+ \
        0.5*gamma_dynamic(gamma, epoch)*tf.reduce_mean(tf.pow(tf.nn.relu(s),2))
        return loss
    
    # Define Gradient    
    def grad(model,x_1,x_2,x_3,epoch):
        with tf.GradientTape() as tape:
            loss_value = loss(model,x_1,x_2,x_3,epoch)
        return loss_value, tape.gradient(loss_value,model.trainable_variables)

    def grad_dynamic(model,x_1,x_2,x_3,p_121,p_131,p_231,epoch):
        with tf.GradientTape() as tape:
            loss_value = loss_dynamic(model,x_1,x_2,x_3,p_121,p_131,p_231,epoch)
        return loss_value, tape.gradient(loss_value,model.trainable_variables)

    # Dynamic Batch Size, Increase over Time
    def batch_size_dynamic(batch_size, epoch):
        return round(((max_iter-epoch)/max_iter)*(batch_size*dynamic_factor)+(epoch/max_iter)*(batch_size))
    
    # Dynamic Gamma, Increase over Time    
    def gamma_dynamic(gamma, epoch):
        return (((max_iter-epoch)/max_iter)*(gamma*dynamic_factor)+(epoch/max_iter)* (gamma))
    
    # Create Optimizer and Model
    optimizer = tf.keras.optimizers.Adam(learning_rate = l_r, beta_1=0.9, beta_2=0.999)
    model = build_model(dynamic_options)
    losses = []

    if dynamic_options:
        # Training Loop
        for epoch in range(int(max_iter)):
            x_1 = next(generate_samples_mu1(batch_size_dynamic(batch_size,epoch)))
            x_2 = next(generate_samples_mu2(batch_size_dynamic(batch_size,epoch)))
            x_3 = next(generate_samples_mu3(batch_size_dynamic(batch_size,epoch)))
            p_121 = tf.random.uniform(shape = [int(batch_size_dynamic(batch_size,epoch)),1])*tf.nn.relu(x_1-upper_epsilon-lower_epsilon-tf.nn.relu(x_1-option_strike1))+lower_epsilon+tf.nn.relu(x_1-option_strike1)
            p_131 = tf.random.uniform(shape = [int(batch_size_dynamic(batch_size,epoch)),1])*tf.nn.relu(x_1-upper_epsilon-lower_epsilon-tf.nn.relu(x_1-option_strike2))+lower_epsilon+tf.nn.relu(x_1-option_strike2)
            p_231 = tf.random.uniform(shape = [int(batch_size_dynamic(batch_size,epoch)),1])*tf.nn.relu(x_2-upper_epsilon-lower_epsilon-tf.nn.relu(x_2-option_strike2))+lower_epsilon+tf.nn.relu(x_2-option_strike2)
            loss_value, grads = grad_dynamic(model, x_1,x_2,x_3,p_121,p_131,p_231, epoch)
            optimizer.apply_gradients(zip(grads, model.trainable_variables))
            losses.append(loss_value.numpy()*(1-2*int(lower_bound)))
            if epoch % 100 == 0 and epoch > 0:
                print("Iteration:{}, Avg. Loss: {}".format((epoch),np.mean(losses[-(round(epoch*0.05))])))
    else:
        # Training Loop
        for epoch in range(int(max_iter)):
            x_1 = next(generate_samples_mu1(batch_size_dynamic(batch_size,epoch)))
            x_2 = next(generate_samples_mu2(batch_size_dynamic(batch_size,epoch)))
            x_3 = next(generate_samples_mu3(batch_size_dynamic(batch_size,epoch)))
            loss_value, grads = grad(model, x_1,x_2,x_3,epoch)
            optimizer.apply_gradients(zip(grads, model.trainable_variables))
            losses.append(loss_value.numpy()*(1-2*int(lower_bound)))
            if epoch % 100 == 0 and epoch > 0:
                print("Iteration:{}, Avg. Loss: {}".format((epoch),np.mean(losses[-(round(epoch*0.05))])))
      
            
    print("Iteration result: {}".format(np.mean(losses[-(round(max_iter*0.05))])))
    return np.mean(losses[-(round(max_iter*0.05))]), model<jupyter_output><empty_output><jupyter_text>Example with lognormal marginals.Compute the grid<jupyter_code>t_1 = 1
t_2 = 2
t_3 = 3
vol = 0.1
S_0 = 1


def log_norm1(batch_size):
    while True:
        x = S_0*tf.math.exp(np.sqrt(vol*t_1)*tf.random.normal(shape = [int(batch_size),1])-vol*(t_1/2))
        yield x
            
def log_norm2(batch_size):
    while True:
        x = S_0*tf.math.exp(np.sqrt(vol*t_2)*tf.random.normal(shape = [int(batch_size),1])-vol*(t_2/2))
        yield x
        
def log_norm3(batch_size):
    while True:
        x = S_0*tf.math.exp(np.sqrt(vol*t_3)*tf.random.normal(shape = [int(batch_size),1])-vol*(t_3/2))
        yield x

def payoff(x,y,z):
   return tf.nn.relu((1/3)*(x+y+z)-S_0)




# Define the grid
epsilon_lower_list  = [e*S_0 for e in [0,0.005,0.01,0.015,0.02]]
epsilon_upper_list = [e*S_0 for e in [0.65,0.7,0.75,0.8,0.85]]

grid = np.zeros([len(epsilon_lower_list),len(epsilon_upper_list)])

for i in range(len(epsilon_lower_list)):
    for j in range(len(epsilon_upper_list)):
        grid[i,j] = compute_bound_3dim(payoff,
                        log_norm1,
                        log_norm2,
                        log_norm3,
                        dynamic_options = True,
                        lower_bound= False,
                          batch_size=2**(13),
                          gamma = 10000,
                          max_magnitude = 1000,
                          max_iter = 50000,
                          nr_neurons = 64*2,
                          depth = 5,
                          l_r =   0.0001,
                          dynamic_factor = 1,
                          option_strike1 = 1,
                          option_strike2 = 1,
                          lower_epsilon = epsilon_lower_list[i],
                          upper_epsilon = epsilon_upper_list[j])[0]
        
grid_df = pd.DataFrame(grid)
grid_df.to_csv("./grid.csv", sep=',',index=False)<jupyter_output><empty_output><jupyter_text>Plotting the result. (Here only for 1 simulation, 30 take a lot of time)<jupyter_code>grid = pd.read_csv("grid.csv") # Here we can also combine more than 1 simulation of grids. In the paper: 30
fig = plt.figure(figsize=(12,8))
ax = Axes3D(fig)

x1 = np.array([0,0.005,0.01,0.015,0.02])*100
#x1 = list(reversed(x1))
y1 = np.array([0.65,0.7,0.75,0.8,0.85])*100
#y1 = list(reversed(y1))
bound_wo = 0.17660294473171234
X, Y = np.meshgrid(x1, y1)
ax.plot_wireframe(X, Y , np.ones((len(x1),len(y1)))*np.max(np.max(average)),color = "red",alpha=1)
ax.plot_surface(X, Y , average,color = "cornflowerblue")
ax.plot_wireframe(X, Y , average,color = "black")
ax.view_init(35,  35)
ax.xaxis.set_rotate_label(False)
ax.yaxis.set_rotate_label(False)
ax.zaxis.set_rotate_label(False)

ax.set_xlabel(r'$\varepsilon_1$', fontsize=25)
ax.set_ylabel(r'$\varepsilon_2$', fontsize=25)
ax.set_zlabel(r'$P_{\Xi_{({p}_{i,j,k}+\varepsilon_1,\overline{p}_{i,j,k}-\varepsilon_2)}}(\Phi)$',
              fontsize=25,rotation = 90)
ax.set_xticks(np.arange(0,2.25,0.25))
ax.set_yticks(np.arange(65,87.5,2.5))
ax.set_zticks([14.5,15,15.5,16,16.5,17])
ax.set_xticklabels([0,0.25,0.5,0.75,1,1.25,1.5,175,2], fontsize=10)
ax.set_yticklabels([65,67.5,70,72.5,75,77.7,80,82.5,85], fontsize=10)
ax.set_zticklabels([14.5,15,15.5,16,16.5,17], fontsize=10)
plt.savefig('3d_improved_bounds.eps', format='eps')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14227752231539087, 'fraction_numerical': 0.06140113605626183, 'mean_word_length': 2.5291123130766784, 'pattern_counts': {'":': 6, '<': 13, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '67243', 'license': 'Unlicense', 'n_tokens_mistral': 4553, 'n_tokens_neox': 4242, 'n_words': 548, 'path': 'Example_3_3.ipynb', 'repo_name': 'juliansester/dynamic_option_trading', 'size': '11091'}