I found a pretty simple solution, which is to use the lemon package.  This cool little extension to ggplot gives you functions facet_rep_grid() and facet_rep_wrap() are the same as facet_grid() and facet_wrap(), respectively, but preserve axis lines and labels on all panels.  This seems pretty much what you want, right OP?
library(lemon)
ggplot(dt, aes(median, sales)) +
  geom_point() +
  facet_rep_grid( month ~ year,
              labeller = labeller( year=label_both, month=label_value),
              repeat.tick.labels = TRUE,
              scales = "free_y") +
  # I recommend increasing panel spacing to account for axis labels
  theme(panel.spacing = unit(20, units="pt"))

As for the coloring each facet according to some scheme... you can always use the trick to draw a geom_rect() behind your data points with some alpha applied.  The key here is that if you use geom_rect() and your original data, you'll get tons of overplotting, so your panels never match the keys in the legend.  The idea is to create a data frame used for geom_rect() that includes only one observation per facet (per combination of year and month in this case.  You also need to include median and sales, since they are designated as x and y in the mapping, but it doesn't really matter the value here.
library(dplyr)
library(tidyr)
ggplot(dt, aes(median, sales)) +
  geom_rect(
    data=expand_grid(month=1:3, year=2000:2002) %>% mutate(sales=0, median=0),
    xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf,
    aes(fill=factor(month)), alpha=0.2
  ) +
  geom_point() +
  facet_rep_grid( month ~ year,
              labeller = labeller( year=label_both, month=label_value),
              repeat.tick.labels = TRUE,
              scales = "free_y") +
  theme(panel.spacing = unit(20, units="pt"))
