<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def plot_frequency_recency_matrix( model, T=1, max_frequency=None, max_recency=None, title=None, xlabel="Customer's Historical Frequency", ylabel="Customer's Recency", **kwargs ):
 """ Plot recency frequecy matrix as heatmap. Plot a figure of expected transactions in T next units of time by a customer's frequency and recency. Parameters model: lifetimes model A fitted lifetimes model. T: fload, optional Next units of time to make predictions for max_frequency: int, optional The maximum frequency to plot. Default is max observed frequency. max_recency: int, optional The maximum recency to plot. This also determines the age of the customer. Default to max observed age. title: str, optional Figure title xlabel: str, optional Figure xlabel ylabel: str, optional Figure ylabel kwargs Passed into the matplotlib.imshow command. Returns ------- axes: matplotlib.AxesSubplot """

    from matplotlib import pyplot as plt

    if max_frequency is None:
        max_frequency = int(model.data["frequency"].max())

    if max_recency is None:
        max_recency = int(model.data["T"].max())

    Z = np.zeros((max_recency + 1, max_frequency + 1))
    for i, recency in enumerate(np.arange(max_recency + 1)):
        for j, frequency in enumerate(np.arange(max_frequency + 1)):
            Z[i, j] = model.conditional_expected_number_of_purchases_up_to_time(T, frequency, recency, max_recency)

    interpolation = kwargs.pop("interpolation", "none")

    ax = plt.subplot(111)
    pcm = ax.imshow(Z, interpolation=interpolation, **kwargs)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    if title is None:
        title = (
            "Expected Number of Future Purchases for {} Unit{} of Time,".format(T, "s"[T == 1 :])
            + "\nby Frequency and Recency of a Customer"
        )
    plt.title(title)

    # turn matrix into square
    forceAspect(ax)

    # plot colorbar beside matrix
    plt.colorbar(pcm, ax=ax)

    return ax