Use numpy. The linear regression code looks something like this:

        A = np.vstack([x, np.ones(len(x))]).T  
        m, c = np.linalg.lstsq(A, y)[0]  

You'll need to import numpy as np.

Here's the documentation: https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.linalg.lstsq.html

Here's how you'd get the x and y in your custom factor:

class linRegres(CustomFactor):  
    inputs = [USEquityPricing.close]  
    def compute(self, today, assets, out, close):  
        x = returns[sid(8554)] #SPY  
        y = pd.DataFrame(close, columns=assets).pct_change()[1:]  
        # ...  

Hope that gets you started. Sorry I can't be of more help. I barely know what I'm doing in python.

did you test this code before posting it as an answer? because it doesn't work

What do you mean it doesn't work? Of course it works. The linalg.lstsq code is straight from the documentation link I posted. Yes, I've tested it. I've been using it.

I'm guessing I have to define "returns" somewhere like

returns = Returns(window_length=2)  

right?

yeah with a customfactor you set something like this:

out[:] = whatever_value_you_want_to_return

not sure how you would return alpha (intercept), beta (slope), r value, p value, and std from one customfactor. I don't know enough python/quantopian to answer that. Somebody else will have to answer.

class linRegres(CustomFactor):  
    outputs = ['alpha', 'beta', 'r_value']  
    inputs = [USEquityPricing.close]  
    def compute(self, today, assets, out, close):  
        x = returns[sid(8554)] #SPY  
        y = pd.DataFrame(close, columns=assets).pct_change()[1:]  
        # ...  

      for i in range(close.shape[1]):  
        out[i].alpha = ...  
        out[i].beta = ...  
        out[i].r_value = ...  

There's also RollingLinearRegressionOfReturns, which might be simpler.