Yes, it is possible to specify 'Upper' as a model that is the upper bound of all terms to consider. For example:
y = 1 + x1 + x2 + x3 + x1.*x2 + x2.*x3 + x1.*x3 + randn(100,1)/10;
d = dataset(x1,x2,x3,y);
GeneralizedLinearModel.stepwise(d,'y~x1+x2+x3')
GeneralizedLinearModel.stepwise(d,'y~x1+x2+x3','Upper','y~x1+x2+x3+x2:x3','verb',2)
The first stepwise invocation will give all pairwise interactions. The second will allow only the specified interaction.
If you have too many predictors to make it feasible to write a formula, you can supply an equivalent terms matrix:
GeneralizedLinearModel.stepwise(d,'y~x1+x2+x3','Upper',[0 0 0 0;1 0 0 0;0 1 0 0;0 0 1 0;0 1 1 0])
This matrix option is described in "help LinearModel.fit".
