Using 'margins, constrast' and marginsplot to test and plot differences across two groups

[Stefano Lombardi]

Suppose you work with the following dataset and perform the regression below with two continuous covaraties (x1 and x2) interacted with a discrete one (year).
You then want to compute predictive margins in correspondence of 2 values of x1 and x2 (the first and last decile of their distribution) and over time.

A minimal working example is the following:

webuse nlswork, clear
keep if year >= 75
rename ttl_exp x1
rename tenure x2

reg ln_wage c.x1##i.year##c.x2

_pctile x1, nq(10)
loc low_x1 = `r(r1)'
loc high_x1 = `r(r9)'
_pctile x2, nq(10)
loc low_x2 = `r(r1)'
loc high_x2 = `r(r9)'

margins, over(year) at(x1=(`low_x1' `high_x1') x2=(`low_x2' `high_x2')) post coeflegend
marginsplot, x(year) ///
    graphregion(color(white)) bgcolor(white) ///
    plot1opts(color(blue*0.5)) ci1opts(color(blue*0.5)) ///
    plot2opts(color(red*0.5)) ci2opts(color(red*0.5)) ///
    plot3opts(color(blue)) ci3opts(color(blue)) ///
    plot4opts(color(red)) ci4opts(color(red)) ///
    legend(order(5 "low x1, low x2" 6 "low x1, high x2" 7 "high x1, low x2" 8 "high x1, high x2") cols(2) region(lstyle(none)))

I would now want to contrast, for each time period, the two red lines values against each other and the two blue lines values against each other, and plot the 2 sets of differences over time (this means, one set of contrasts "low x1, high x2" vs. "high x1, high x2" over time and another set of contrasts "low x1, low x2" vs. "high x1, low x2" over time).

Then I would like to have a separate set of contrasts of the type "mild red vs. mild blue" and "deep red vs. deep blue" over time in a separate graph.

I have looked into "margins, contrast" but with no success. Any help is much appreciated.