Interaction terms in regression when variables can be negative Suppose I have a regression
$y = b_0 + b_1 x_1 + b_2 x_2 $
where both $x_1$ and $x_2$ have a range from $-\infty$ to $\infty$ and have been centered.
The correlation between $x_1$ and $x_2$ is $0.3$.
Visually, it seems like both $x_1$ and $x_2$ are roughly almost linear in $y$ except in the extremes, with a lot of noise.
I want to add an interaction variable $x_1 * x_2$, but the issue is when both $x_1$ and $x_2$ are negative, their multiplication becomes positive, which is not the desired effect, as the effect I want to capture is that when both are strongly positive, $y$ is more positive than what the linear sum would give, and similarly when both are strongly negative, $y$ is more negative than what the linear sum would give.
I googled online and there doesn’t seem to be any good reference on how to handle (create) interaction variables when both have a real number domain.
Is there any good reference or recommendations to create interaction variables when both $x_1$ and $x_2$ have a range from $-\infty$ to $\infty$ ?
 A: Maybe you can think of using a dummy to pre-multiply the interaction term so that the effect of the interaction term on y is not considered for cases where x1 and x2 are both negative (as you set the dummy to go to 0 in that case). Or it is only considered for cases where both x1 and x2 are strictly positive (which seems to be the case of your interest, right? It’s up to you). Let’s denote the former case as A and the latter as B.
So you would have an interaction term like $\beta \gamma x_{1} x_{2} $ where the dummy is $\gamma$. Under case A, $\gamma$ is 0 if x1 and x2 are both negative and 1 otherwise (so you just dismiss the troubling cases where both the regressors are negative). Instead, under case B, if you just want to assess the effect of cases where both x1 and x2 are positive and dismiss what happens in all other circumstances (including the mixed cases where one between x1 and x2 is negative), then $\gamma$ must be set to 1 iff x1>0 and x2>0, and 0 otherwise.
A: I think, your problem is not the interaction term per se but that you are actually not interested in an interaction model but in a model of a different shape. You would need to tell us more about your actual hypothesis to get helpful concrete suggestions but my impression is that your question could be analyzed in the framework of response surface analysis which is a regression model that includes more non-linear terms beyond the "simple" interaction. It is especially suited to analyze the consequences of certain combinations of two predictor variables on the dependent variable.
Here's an R package for that:
https://github.com/nicebread/fSRM
You can also find many references there.
A: This is not about the range of the variables in the interaction, but that you seem to want a nonlinear interaction. That is better framed as a nonparametric model of the regression function, which could be represented via splines.  I will simulate some data with some nonlinear interaction, and use R's mgcv package to fit a regression represented by a thin plate spline, see for instance Smoothing methods for gam in mgcv package?.  
The fitted regression function can be represented with a contour plot:

where also uncertainty bands is shown around the contours. 
The R code used is:
set.seed(7*11*13) # my public seed 
N  <- 1000
x1 <- rnorm(N, 0, 10)
x2 <- rnorm(N, 0, 10)
inter <- ifelse( (x1<0)&(x2<0), -1*0.05*x1*x2, 0.05*x1*x2 )
mu <- 0 + 0.1*x1 + 0.1*x2 + inter
y  <- mu + rnorm(N, 0, 5)  

mydf <- data.frame(x1, x2, mu, inter, y)       

library(splines)
library(mgcv)  

mod.gam <- mgcv::gam(y ~ s(x1, x2, bs="tp"), data=mydf)

summary(mod.gam)
Family: gaussian 
Link function: identity 

Formula:
y ~ s(x1, x2, bs = "tp")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -2.0795     0.1582  -13.14   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df     F p-value    
s(x1,x2) 25.23  28.18 49.24  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.582   Deviance explained = 59.2%
GCV = 25.708  Scale est. = 25.034    n = 1000

plot(mod.gam)

For comparison, some simpler models could be fit with
mod0 <- lm(y ~ x1 + x2, data=mydf)
mod1 <- lm(y ~ x1 + x2 + x1:x2, data=mydf)

mod2 <- lm(y ~ x1 + x2 + I(x1^2) + I(x2^2) + x1:x2, data=mydf)

mod.spline <- lm(y ~ x1 + x2 + ns(x1 * x2, df = 6), data=mydf)    

