How to calculate joint feature contribution for XGBoost Classifier in python?
I did some research and learnt about xgbfir package. It gives the joint contributions into an excel file. You can set the level of interaction with this. I wrote some code around it to generate a plot that solves the purpose. 

If the package is not installed

pip install xgbfir


After the installation:

import xgbfir
from matplotlib import pyplot as plt

xgbfir.saveXgbFI(model, feature_names=X.columns, OutputXlsxFile='FI.xlsx')

joint_contrib = pd.read_excel('FI.xlsx')

xls = pd.ExcelFile('FI.xlsx')
df1 = pd.read_excel(xls, 'Interaction Depth 0')
df2 = pd.read_excel(xls, 'Interaction Depth 1')
df3 = pd.read_excel(xls, 'Interaction Depth 2')

frames = [df1, df2, df3]
joint_contrib = pd.concat(frames)

joint_contrib=joint_contrib.sort_values(by='Gain', ascending=True)
joint_contrib=joint_contrib.head(20)

height = joint_contrib['Gain']
bars = joint_contrib['Interaction']
y_pos = np.arange(len(bars))

plt.barh(y_pos, height)
plt.yticks(y_pos, bars)
plt.show()


This will give the top 20 feature interactions in terms of gain. 

Thanks to Philip Cho who introduced me to xgbfir. 

Follow the link for more information regarding xgbfir