Nonlinear Autoregressive model parameter estimation from time series I'm working on a nonlinear multivariate autoregressive model of order 1 (markovian).
It is a discrete-time dynamical system which models exchange of mass between compartments in a compartmental model where each compartment corresponds a given population type characterized, for every time instant, by a fraction of mass $x_u(t)$. The network dynamics are expressed by quantities of mass which have chances of interact and being converted to another type, thus leading to coupled interactions which determine flow of mass between types.
The mass for type $u$ at time $t + 1$ is given by its mass at $t$ plus the incoming flows, coming from a subset of neighboring types, minus outgoing flow to other types.
$$
x_u(t+1) = x_u(t) + \sum_{v : u \in N_{\mathcal C}(v)} f_{vu}(t) -\sum_{v \in \mathcal N_{\mathcal C}(u)} f_{uv}(t)
$$
Where each flow has the form:
$$
f_{uv}(t) = x_u (t) \cdot \sum_{z \in S_u} \frac{x_z(t)}{(\alpha_{u} - 1 )\cdot x_u(t) + B
}  \cdot \frac{p_{uz}}{m}
$$
My objective is doing estimation of parameters $\alpha_u$ (1 for each compartment) and $p_{uv}$ (one for each couple of compartments) from a time series of observations concerning the masses $x_u$. I obtained good results using least squares nonlinear parameters estimation for low dimensional systems, while it quickly becomes impractical for systems with moderate to high number of compartments. Could someone suggest a technique/direction to learn the parameters of this kind of models?
EDIT: the parameters $p_{uv}$ are constrained to be probabilities (thus ranging in [0,1]) and $\alpha_u \geq 1$ 
 A: Let's denote $\hat{x}_u(t)$ the mass predicted by the model given the parameter estimates $\hat{\alpha}_u$ and $\hat{p}_{uv}$. Then, you can define the following function:
$$
S(\alpha,p_{uv};x_u) = \sum_{u,t} \left( x_u(t) - \hat{x}_u(t) \right)^2 \,.
$$
This function returns a scalar that is the sum over all $u$ and $t$ of the squared differences between the observed values and the predicted values. 
Thus, it is a measure of the prediction errors given a set of parameter estimates. (The square is taken to avoid that negative errors are offset positive errors.)
Then, the goal is to minimize this error function. This could be done by means of a general purpose optimization algorithm, for example: the Nelder-Mead algorithm
or the Broyden–Fletcher–Goldfarb–Shanno algorithm.
Starting from  an arbitrary set of parameters, these algorithms follow an iterative procedure seeking for the values of the parameters that minimises the objective function. The values of the parameters found by this approach will be optimal in the sense that minimise the prediction errors.
Note: If there is more than one local minima, there is no guarantee that these algorithms will converge to a global minimum. Running the algorithms from different starting values is advisable.
