好的，你首先要做的就是绘制数据。为了简单起见，我复制了this figure，因为它只有8个事件发生，所以很容易看到系统的行为。以下代码：import numpy as np

import math, matplotlib

import matplotlib.pyplot

import matplotlib.lines

mu = 0.1 # Parameter values as found in the article http://jheusser.github.io/2013/09/08/hawkes.html Hawkes Process section.

alpha = 1.0

beta = 0.5

EventTimes = np.array([0.7, 1.2, 2.0, 3.8, 7.1, 8.2, 8.9, 9.0])

" Compute conditional intensities for all times using the Hawkes process. "

timesOfInterest = np.linspace(0.0, 10.0, 100) # Times where the intensity will be sampled.

conditionalIntensities = [] # Conditional intensity for every epoch of interest.

for t in timesOfInterest:

conditionalIntensities.append( mu + np.array( [alpha*math.exp(-beta*(t-ti)) if t > ti else 0.0 for ti in EventTimes] ).sum() ) # Find the contributions of all preceding events to the overall chance of another one occurring. All events that occur after t have no contribution.

" Plot the conditional intensity time history. "

fig = matplotlib.pyplot.figure()

ax = fig.gca()

labelsFontSize = 16

ticksFontSize = 14

fig.suptitle(r"$Conditional\ intensity\ VS\ time$", fontsize=20)

ax.grid(True)

ax.set_xlabel(r'$Time$',fontsize=labelsFontSize)

ax.set_ylabel(r'$\lambda$',fontsize=labelsFontSize)

matplotlib.rc('xtick', labelsize=ticksFontSize)

matplotlib.rc('ytick', labelsize=ticksFontSize)

eventsScatter = ax.scatter(EventTimes,np.ones(len(EventTimes))) # Just to indicate where the events took place.

ax.plot(timesOfInterest, conditionalIntensities, color='red', linestyle='solid', marker=None, markerfacecolor='blue', markersize=12)

fittedPlot = matplotlib.lines.Line2D([],[],color='red', linestyle='solid', marker=None, markerfacecolor='blue', markersize=12)

fig.legend([fittedPlot, eventsScatter], [r'$Conditional\ intensity\ computed\ from\ events$', r'$Events$'])

matplotlib.pyplot.show()

即使我选择的事件时间有点武断：

这也可以应用于5000笔交易的example set of data集合，方法是将数据绑定并将每个bin视为一个事件。然而，现在的情况是，由于每个箱子中发生的交易数量不同，每个事件的权重略有不同。

这也在the article将比特币交易到达拟合到霍克斯过程中的一节中提到，并提出了解决这个问题的方法：The only difference to the original dataset is that I added a random millisecond timestamp to all trades that share a timestamp with another trade. This is required as the model requires to distinguish every trade (i.e. every trade must have a unique timestamp).这一点包含在以下代码中：

^{pr2}$

这将生成以下拟合图：

所有这些看起来都不错，但是，当你仔细观察细节时，你会发现，仅仅通过取交易数量的一个向量并减去拟合的向量来计算残差是行不通的，因为它们的长度不同：

然而，有可能提取与经验数据记录相同时期的强度，然后计算残差。这使您能够找到经验数据和拟合数据的分位数，并将它们相互对比，从而生成QQ图：""" GENERATE THE QQ PLOT. """

" Process the data and compute the quantiles. "

orderStatistics=[]; orderStatistics2=[];

for i in range( empirical_1min.values.size ): # Make sure all the NANs are filtered out and both arrays have the same size.

if not np.isnan( empirical_1min.values[i] ):

orderStatistics.append(empirical_1min.values[i])

orderStatistics2.append(conditionalIntensities[i])

orderStatistics = np.array(orderStatistics); orderStatistics2 = np.array(orderStatistics2);

orderStatistics.sort(axis=0) # Need to sort data in ascending order to make a QQ plot. orderStatistics is a column vector.

orderStatistics2.sort()

smapleQuantiles=np.zeros( orderStatistics.size ) # Quantiles of the empirical data.

smapleQuantiles2=np.zeros( orderStatistics2.size ) # Quantiles of the data fitted using the Hawkes process.

for i in range( orderStatistics.size ):

temp = int( 100*(i-0.5)/float(smapleQuantiles.size) ) # (i-0.5)/float(smapleQuantiles.size) th quantile. COnvert to % as expected by the numpy function.

if temp<0.0:

temp=0.0 # Avoid having -ve percentiles.

smapleQuantiles[i] = np.percentile(orderStatistics, temp)

smapleQuantiles2[i] = np.percentile(orderStatistics2, temp)

" Make the quantile plot of empirical data first. "

fig2 = matplotlib.pyplot.figure()

ax2 = fig2.gca(aspect="equal")

fig2.suptitle(r"$Quantile\ plot$", fontsize=20)

ax2.grid(True)

ax2.set_xlabel(r'$Sample\ fraction\ (\%)$',fontsize=labelsFontSize)

ax2.set_ylabel(r'$Observations$',fontsize=labelsFontSize)

matplotlib.rc('xtick', labelsize=ticksFontSize)

matplotlib.rc('ytick', labelsize=ticksFontSize)

distScatter = ax2.scatter(smapleQuantiles, orderStatistics, c='blue', marker='o') # If these are close to the straight line with slope line these points come from a normal distribution.

ax2.plot(smapleQuantiles, smapleQuantiles, color='red', linestyle='solid', marker=None, markerfacecolor='red', markersize=12)

normalDistPlot = matplotlib.lines.Line2D([],[],color='red', linestyle='solid', marker=None, markerfacecolor='red', markersize=12)

fig2.legend([normalDistPlot, distScatter], [r'$Normal\ distribution$', r'$Empirical\ data$'])

matplotlib.pyplot.show()

" Make a QQ plot. "

fig3 = matplotlib.pyplot.figure()

ax3 = fig3.gca(aspect="equal")

fig3.suptitle(r"$Quantile\ -\ Quantile\ plot$", fontsize=20)

ax3.grid(True)

ax3.set_xlabel(r'$Empirical\ data$',fontsize=labelsFontSize)

ax3.set_ylabel(r'$Data\ fitted\ with\ Hawkes\ distribution$',fontsize=labelsFontSize)

matplotlib.rc('xtick', labelsize=ticksFontSize)

matplotlib.rc('ytick', labelsize=ticksFontSize)

distributionScatter = ax3.scatter(smapleQuantiles, smapleQuantiles2, c='blue', marker='x') # If these are close to the straight line with slope line these points come from a normal distribution.

ax3.plot(smapleQuantiles, smapleQuantiles, color='red', linestyle='solid', marker=None, markerfacecolor='red', markersize=12)

normalDistPlot2 = matplotlib.lines.Line2D([],[],color='red', linestyle='solid', marker=None, markerfacecolor='red', markersize=12)

fig3.legend([normalDistPlot2, distributionScatter], [r'$Normal\ distribution$', r'$Comparison\ of\ datasets$'])

matplotlib.pyplot.show()

这将生成以下绘图：

^{5美元

经验数据的分位数图不完全是the same as in the article，我不知道为什么，因为我不擅长统计学。但是，从编程的角度来看，这就是您可以完成所有这些的方法。在