Estimation of a system Suppose we have a system that essentially evolves as follows:
stock_t+1  = stock_t + inflows_t - outflows_t 
inflows_t  = a1*predictor11_t + a2*predictor12_t+.... error1_t
outflows_t = b1*predictor21_t + b2*predictor22_t+.... error2_t

I have observations for each of these variables, i.e. I can observe the stock, the inflows, outflows, and the sum of the flows (they add up correctly), as well as each of the predictors. All variables are time series, and simple time series analysis goes a long way. That said, while the stock is not overly volatile and relatively easy to predict, the flow variables are much more volatile and more challenging to model.
Using regression analysis, in the beginning I've only attempted to model the evolution of the stock using a combination of the flow predictors. I found that the best predictors of the stock are close to the best predictors of the flows (though not identical, i.e. some additional transformation is required). Recently, I have also attempted to model the flows. 
The trouble is - not unexpectedly, I should say - that using the best regression models for each of the series, the system does not follow the add-up constraint; or put differently: if I calculate the evolution of the stock using the first equation, based on a starting value and equations 2-3, then my prediction of the stock is quite a bit different than my direct forecast of the stock. At the same time, since the stock model is the model that I have the most faith in, I'd rather not move away from these predictions.
So I was thinking that there surely must be a way to model the entire system directly, rather than estimate each equation separately. Right now a state-space approach comes to mind; before I go off into that direction, though, I am wondering whether I am missing something and whether anyone has a different suggestion.
PS I'm using R
 A: First let's give some familiar structure and definitions to your formulation:
$s_{t+1} = s_{t} + \Delta_t$
$x_t = \mathbf{a}^\intercal \mathbf{p}_t^{(1)} + \epsilon^1_t$
$y_t = \mathbf{b}^\intercal \mathbf{p}_t^{(2)} + \epsilon^2_t$
$\Delta_t = x_t - y_t$. 
So the "flows" $x$ and $y$ have the structure of inner products with additive noise (what distribution?)
I'm not familiar with the assumptions made in economics, nor do they matter here for now. (you mention "volatility", which is meaningless in mathematics, until you attach distributions and parameters to all your stochastic processes)
So you say you have to predict $s_{t+1}$ based on $x_t$ and $y_t$? What about the other parameters, e.g. the noise terms $\epsilon$, the coefficients $\mathbf{a}$ and $\mathbf{b}$ and the (time-varying) "predictors" (e.g. basis functions) $\mathbf{p}_t^{(i)}$?
I'm afraid the problem specification, as you gave it, is incomplete.
If you observe discrepancies between two approaches, it's best if you start by showing us all the models you use, in a unified notation.
