What is the reason for I(1) integration order limit of independent (or dependent) variable in ARDL regression?, to be specific I(2) variable will 'break' the ARDL model/estimation. Fast thinking I cannot see where the issue is in estimation - e.g. how the model breaks, my guess is that the model's residuals would be then non-stationary?
Also what model you would suggest to use if one of time-series being modeled is of I(2) order or higher? Dynamic Ordinary Least Squares Estimator?
...I(1) integration order limit of independent (or dependent) variable in ARDL regression...
It's not specified exactly in what context this statement is made, but in principle there's no such "limit".
Inference can be carried out for any model for which the asymptotic distributions can be obtained. In the case of ARDL model with I(2) or higher variables, the asymptotic distributions are non-normal but still accessible.
Consider the simplest such model $$ w_t = w_{t-1} + \epsilon_t, \;\; \epsilon_t = \epsilon_{t-1} + u_t, \;\; u_t \stackrel{i.i.d.}{\sim} (0, \sigma^2). $$ There is an autoregression of an I(2) variable $w_t$ and I(1) error term $\epsilon_t$.
It is clear, at least heuristically, that the following is true:
$$
\frac{1}{T^{\frac12} \sigma} \sum_{t=1}^{[rT]} \epsilon_t \Rightarrow W(r), \;\;
\mbox{and}\;\;
\frac{1}{T^{\frac{3}{2}} \sigma}\sum_{t=1}^{[rT]} w_t \Rightarrow \int_0^r W(s) ds = X(r)
$$
where $W$ is a standard Brownian motion and "$\Rightarrow$" denotes weak convergence of stochastic processes on the unit interval [0,1].
For the OLS estimate $\hat{\phi}$ from regressing $w_t$ on $w_{t-1}$, we then have $$ T(\hat{\phi} - 1) = T \cdot\frac{\sum_{t = 1}^T \epsilon_t w_{t-1} }{\sum_{t = 1}^T w_{t-1}^2} =\frac{ \frac{1}{T^3}(w_T^2 - w_0^2 - \frac12 \sum_{t=1}^T \epsilon_t^2) }{\frac{1}{T^4}\sum_{t = 1}^T w_{t-1}^2} \Rightarrow \frac{X^2(1)}{\int_0^1X (r) dr}. $$
So $\hat{\phi}$ is (super-)consistent, converging to true parameter 1 at rate $\frac{1}{T}$. This would remain true under very general dependent structure for $u_t$. The i.i.d. assumption is not necessary.
As a general rule, OLS estimate is (super-)consistent, and $R^2$ approaches 1 as sample size gets large, whenever the regressors have higher order of integration than the error term. You can also observe this in, for example, a cointegration regression.
...my guess is that the model's residuals would be then non-stationary?
Yes, the residuals would be non-stationary, since the error term $\epsilon_t$ is non-stationary, but this is again not necessarily a problem.
(The residual sum of squares from the regression is $$ \sum_{t=1}^T \epsilon_t^2 - \frac{\left( \sum_{t=1}^T \epsilon_t w_{t-1} \right)^2}{\sum_{t=1}^T w_{t-1}^2} = O_p(T^2) - O_p(T^2) = O_p(T^2). $$ This tells you that the residual is I(1).)
By similar kind of algebra, one can obtain the asymptotic distribution for the usual t-statistic (with scaling by $\sqrt{T}$) $\frac{1}{\sqrt{T}} \frac{\hat{\phi} - 1}{\mbox{s.e.}}$, where s.e. is the standard error.
what model...to use if one of time-series being modeled is of I(2) order or higher?
Above discussion shows that, in principle, you can carry out inference for ARDL models with variables with integration order higher than 1.
Whether you want to actually fit such models in practice is a different question. The Consumer Price Index is one common series that is, arguably, I(2), but I have never seen it fitted as such. Given that there are far more known techniques for dealing with I(0) and I(1) series than series of higher integration order, I would suggest differencing to lower the integration order if possible.
