How do I test whether two regressions lines cross at y = 0? I have measurements for daily space heating vs daily mean outdoor temperature for two different control strategies. The data is shown here:

I have also performed linear regression, of the form:
lm(Energy ~ Control * MeanOutdoorTemp)

This yields four coefficients:
Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                170.8293     4.2083  40.594  < 2e-16 ***
ControlMPC                 -30.7044     4.9025  -6.263 1.38e-07 ***
MeanOutdoorTemp             -6.2924     1.6466  -3.821 0.000413 ***
ControlMPC:MeanOutdoorTemp   0.8211     1.7162   0.478 0.634709    

I'd like to test whether these two regression lines cross at zero energy. What would be the statistically correct way to do that?
 A: The equation for the 1st line is $\beta_0 +\beta_1*x$ and part of the question is where this value equals 0, so you could set this equal to 0 and solve for $x$ to get $x=\frac{-\beta_0}{\beta_1}$.  The equation for the 2nd line would then be $(\beta_0 + \beta_2) + (\beta_1 + \beta_3) * x$ and you can again set this to 0 and solve for $x$ to get $x=\frac{-(\beta_0 + \beta_2)}{(\beta_1 + \beta_3)}$.  You can now test if the difference between these 2 is 0, or fit a confidence interval to the difference.
Unfortunately, working with ratios in cases like this can be very difficult (the distribution of a ratio of 2 normals can be a Cauchy if the denominator can be close to 0).  One approach to get a reasonable approximation is to use simulation.  This will assume that the estimates are normally distributed (standard assumption, but make sure that your residuals are normal enough or sample size large enough for this to be reasonable)  You will need the covariance of your slope estimates (vcov(fit) in R where fit is the fitted model object).  Now you can simulate slopes from a multivariate normal with mean equal to the estimated slopes and the given covariance matrix (in R use the mvrnorm function from the MASS package).  The for each simulated set of values (4 slopes) compute the difference of the 2 equations above.  Do this for a few thousand simulated values and then plot/summarize the results and see how they compare to 0.
Here is some sample R code to do a simulation:
x <- runif(50, 5, 10)
g <- rep(0:1, each=25)
y <- 0.2*x + 0.1*x*g + rnorm(50, 0, 0.1)
# y <- 0.2*x - 0.5*g + 0.1*x*g + rnorm(50, 0, 0.1)
plot(x,y, pch= g+1, xlim=c(0,10), ylim=range(0,y))

fit <- lm(y~x*g)
summary(fit)

tmp <- coef(fit)
abline( tmp[1], tmp[2], col='blue' )
abline( tmp[1] + tmp[3], tmp[2] + tmp[4], col='green' )

library(MASS)
sims <- mvrnorm(100000, tmp, vcov(fit))
diffs <- sims[,1]/sims[,2] - (sims[,1]+sims[,3])/(sims[,2]+sims[,4])
hist(diffs)
abline(v=0, col='red')
mean(diffs>0)
mean(diffs<0)
quantile(diffs, c(0.025, 0.975))

You can play around with this (the commented line gives different x intercepts by changing the y intercept of the 2nd group, run this line instead to see a difference, the uncommented line they do have the same x intercept (at 0)).
Try playing around with the parameters to see the effects.  This can also give you a feel for power and what size of difference you can detect.
