For instance, I've got K variables and want to select the optimal lag structure by searching across the (pmax-1)^K ARDL models, spanned by p=0,1,..., pmax and p_k=0,1,..., pmax, k=1,.., K using some criteria (e.g. AIC).

I am looking for a clever trick to estimate all of these possible models but only have a complicated "loop-in-loop-in-loop" structure in mind. Does anybody have a clever idea how to realize this task efficiently?