# Maximum Entropy Bootstrap - how to implement?

I have a question about the Maximum Entropy Bootstrap algorithm for time series.

I'm confused about how to implement step 5, where you compute the sample quantiles of the Maximum Entropy distribution. Essentially, how do you go from the red stars, to the blue triangles on the left in this diagram:

The steps of the algorithm are shown below:

1. Sort the original data in increasing order to create order statistics x(t) and store the ordering index vector.
2. Compute intermediate points zt = (x(t) + x(t+1))/2 for t = 1, . . . , T − 1 from the order statistics. Hrishikesh D. Vinod, Javier L´opez-de-Lacalle 3
3. Compute the trimmed mean mtrm of deviations xt − xt−1 among all consecutive observations. Compute the lower limit for left tail as z0 = x(1) − mtrm and upper limit for right tail as zT = x(T) + mtrm. These limits become the limiting intermediate points.
4. Compute the mean of the maximum entropy density within each interval such that the ‘mean-preserving constraint’ (designed to eventually satisfy the ergodic theorem) is satisfied. Interval means are denoted as mt . The means for the first and the last interval have simpler formulas.
5. Generate random numbers from the [0, 1] uniform interval, compute sample quantiles of the ME density at those points and sort them.
6. Reorder the sorted sample quantiles by using the ordering index of step 1. This recovers the time dependence relationships of the originally observed data.
7. Repeat steps 2 to 6 several times (e.g., 999).

The diagram comes from: https://cran.r-project.org/web/packages/meboot/vignettes/meboot.pdf

I've written some python code to do these steps, but I'm not sure it's correct:

def meboot(x, col):
xt = x.sort_values(by=[col])
zt = xt.rolling(2).mean()
diff = x.rolling(window=2).apply(lambda x: x[1] - x[0])
trimmed_mean = diff.mean()
upper = xt.iloc[-1] + trimmed_mean
lower = xt.iloc[0] - trimmed_mean
zt.iloc[0] = lower
zt.iloc[-1] = upper

desired_means = pd.DataFrame(data=np.zeros(x.shape[0],))

desired_means.iloc[0] = (0.75*xt.iloc[0] + 0.25*xt.iloc[1]).values

for k in range(1, xt.shape[0]-1):
desired_means.iloc[k] = (0.25*xt.iloc[k-1] + 0.5*xt.iloc[k] + 0.25*xt.iloc[k+1]).values

desired_means.iloc[-1] = (0.75*xt.iloc[-1] + 0.25*xt.iloc[-2]).values

U = pd.DataFrame(np.sort(np.random.rand(xt.shape[0])), columns=[col], index=xt.index)

N = xt.shape[0]
x = [float(y)/N for y in range(N+1)]
inds = range(len(x))
quantiles = np.zeros(N)

for k in range(N):
yy = np.abs(x - U.iloc[k].values)
ind = np.argmin(yy)
if x[ind] > U.iloc[k].values:
ind -= 1
c = (2*desired_means.iloc[ind-1].values - zt.iloc[ind-1].values - zt.iloc[ind].values) / 2
y0 = zt.iloc[ind - 1].values + c
y1 = zt.iloc[ind].values + c
quantiles[k] = y0 + (U.iloc[k].values - x[ind]) * (y1 - y0) / (x[ind + 1] - x[ind])
quantiles = pd.DataFrame(quantiles, columns=['Passengers'], index=xt.index)
return quantiles.sort_index()

• You may send an email to H. D. Vinod asking for a document he wrote clarifying this toy example. – javlacalle Jul 8 '18 at 20:55
• Hi, have you managed to make it work ? I would be interested to see the final code – jmague Feb 5 at 13:07