Question: Please forgive my ignorance. All I'm trying to do is add a squared term to my regression without going through the trouble of defining a new column in my dataframe. I'm using statsmodels.formula.api (as stats) because the format is similar to R, which I am more familiar with. 

```
hours_model = stats.ols(formula='act_hours ~ h_hours + C(month) + trend', data = df).fit()

```

The above works as expected.

```
hours_model = stats.ols(formula='act_hours ~ h_hours + h_hours**2 + C(month) + trend', data = df).fit()

```

This omits h\_hours\*\*2 and returns the same output as the line above. 

I've also tried: h\_hours^2, math.pow(h\_hours,2), and poly(h\_hours,2) 
All throw errors. 

Any help would be appreciated.

Answer: 
You can try using `I()` like in R:

```
import statsmodels.formula.api as smf

np.random.seed(0)

df = pd.DataFrame({'act_hours':np.random.uniform(1,4,100),'h_hours':np.random.uniform(1,4,100),
                  'month':np.random.randint(0,3,100),'trend':np.random.uniform(0,2,100)})

model = 'act_hours ~ h_hours + I(h_hours**2)'
hours_model = smf.ols(formula = model, data = df)

hours_model.exog[:5,]

array([[ 1.        ,  3.03344961,  9.20181654],
       [ 1.        ,  1.81002392,  3.27618659],
       [ 1.        ,  3.20558207, 10.27575638],
       [ 1.        ,  3.88656564, 15.10539244],
       [ 1.        ,  1.74625943,  3.049422  ]])

```
Currently, although the statsmodels formula API (in fact Patsy library) doesn't support `poly(variable, degree)` function as in R, NumPy's `vander(variable, degree+1)` can do the job. However, pay **attention** that `np.vander()` produces the Vandermonde matrix which means you get intercept column too! Let's see this function in an example:

```
>> x = np.array([1, 2, 3, 5])
>> np.vander(x, 4, increasing=True)

array([[  1,   1,   1,   1],
       [  1,   2,   4,   8],
       [  1,   3,   9,  27],
       [  1,   5,  25, 125]])

```

So, you need to remove Patsy's internal intercept by adding `-1` to your formula:

```
hours_model = stats.ols(formula='act_hours ~ np.vander(h_hours, 3, increasing=True) - 1', data = df).fit()

```

Note that you need to pass `your_desired_degree + 1` because the first column is x^0=1.