import numpy as np
res = 100
T = 1
dt = T / res
kappa_int = dt * np.arange(res)  # [int_0^t kappa(s)ds, for all t]
X = np.ones(res)  # [X[t], for all t]
integrand1 = np.zeros(res)
integrand2 = np.zeros(res)
for u in range(0, res):
    integrand1[u] = np.exp(-kappa_int[u])
    integrand3 = np.zeros(u)
    integrand4 = np.zeros(u)
    for tau in range(0, u):
        integrand3[tau] = np.exp(-kappa_int[tau])
        integrand5 = np.zeros(tau)
        for s in range(0, tau):
            integrand5[s] = np.exp(2*kappa_int[s])*X[s]
        integrand4[tau] = np.sum(integrand5) * dt
    integrand2[u] = np.dot(integrand3, integrand4) * dt
output = 6 * np.dot(integrand1, integrand2) * dt