I am currently working on a class project, in which I need to identify groups of regions in my sample base on a regression tree. However, I am unfortunately not very familiar with regression trees in general, which is why I am a bit confused about their implementation in statistical software packages, such as R or MATLAB.

The algorithm I intend to use, however, is pretty straightforward and is taken from P.Postiglione et al. (2010) - A regression tree algorithm for the identification of convergence clubs. Maybe someone here more, who knows more about regression trees than me can give me advice as to which functions in R or MATLAB might be useful for achieving this, and on whether I will need to implement the algorithm from scratch by myself.

The algorithm is the following (Postiglioni et al. 2010:2778):

1. A club equivalent to the entire population of European regions, say $P = \{i_{1}, i_{2}, {...}, i_n\}$ is generated and the model under investigation estimated. This can be referred to as step $0$ of our procedure. Let us consider the set of clubs generated at step $k \geq 0,$ say $S_k$. For each club $C$ in $S_k$ the following procedure is executed.

2. Let $X(i_j)$ be the value assumed by a given splitting variable $X$ at region $i_j$, and let $X_A$ the set of values of $X$ observed in $A$, i.e. $X_A = \{X(i_j) : i_j \in A\}$. If X is an ordered splitting variable, for every $x \in X_A \backslash \{max(X_a)\}$ the current club $C$ is bi-partitioned in sets $B(x) = \{i_j \in C : X(i_j) \leq x\}$ and $\bar{B}(x) = C \backslash B(x)$. If $X$ is an unordered splitting variable, for every proper subset $B$ of $C$, $C$ is instead bi-partitioned in sets $B$ and $\bar{B} = C\backslash B$.

3. The model under study is thus estimated making use of the sampling information corresponding to both $B$ and $\bar{B}$, say obtaining the Maximum Likelihood estimates $\mathbf{\theta_B}$ and $\mathbf{\theta_\bar{B}}$, respectively. The distance among such parameter vectors is evaluated and inspected for statistical significance using the corresponding probability value. The statistical test is baed upon the following arguments. Provided that $\theta_B$ and $\theta_\bar{B}$ are independent, given their asymptotical normal distributions, we have that the statistic: $$ S = (\theta_B - \theta_\bar{B})^T(\Sigma_B + \Sigma_\bar{B})^{-1}(\theta_B - \theta_\bar{B})$$ follows a Chi Squared distribution with $d$ degrees of freedom, $d$ being the size of $\theta$. The statistic [...] is the objective function of our regression tree.

So, as I mentioned, the intuition behind the algorithm is pretty straight forward. The splitting is then continued until either (i) the last probability value exceeds a given benchmark probability $\tilde{p}$, (ii) the sub-clubs get too small, or (iii) we have reached a desired maximum number of clubs.

I am sorry for this basic question, but honestly I got a bit confused about this whole field of regression trees, although the intuition behind this particular case was pretty clear to me. I just didn't know whether I could specify an objective function like this for the existing tree packages. Also, if I have to implement this one myself, I would be thankful for any functions I could exploit in R or MATLAB to make life easier for myself.

Thanks in advance!

EDIT: I think coding this thing on my own for the case of unordered variables would also be a lot easier if I had a method to a) get all possible combinations of observations if splitting my set of observations $C$ into two distinct subsets $B$ and $\bar{B} = C\backslash B$, since then I could simply store the regression results and use them to compute the statistic above and b) keep track of the observations during the whole process, but my search so far hasn't really yielded anything that I found useful. So already finding help for these sub-problems would be really great.