Predicting Y from a regression model for dY I have some time series data where I'm modelling temperature as a function of various predictors. On physical grounds, I can expect that
$$\frac{dT}{dt} \propto T_a - T$$
where $T_a$ is the ambient temperature (which can vary over time, but whose values are known). I thus fit models of the form
$$\Delta T(t) \sim \alpha + \beta \left[ T_a(t) -T(t) \right] + \gamma X(t)$$
with $X$ being the other covariates, and $\alpha$, $\beta$ and $\gamma$ are the regression parameters. I can fit these easily enough in R:
lm(diff(T) ~ I(Ta - T) + x, data=df)

and I can get predictions for the change in $T$. However, what I really want are predictions for $T$ itself. At the moment I'm calculating these via a loop, where I plug $\hat{T}(t)$ into the regression equation to obtain $\hat{\Delta T}(t+1)$.
Is there any R package, probably time series-related, that will do these calculations automatically?
Also, if there are any issues with this approach, I'd be happy to know about them.
 A: Your approach can be related to using a time-series model either in
discrete time (DT) or continuous time (CT). It seems better to work
with the difference $Z(t) := T(t) - T_{\textrm{a}}(t)$ because
$T_{\textrm{a}}(t)$ will be required as an offset term for forecasting
as well as for fitting. A possible model writes in CT as
$$
\frac{\textrm{d}}{\textrm{d}t}Z(t) + \theta Z(t) = \alpha +
\boldsymbol{\gamma}^\top \mathbf{X}(t) + \varepsilon(t)
$$
where $\varepsilon(t)$ stands for a white noise. The constraint
$\theta > 0$ should hold to ensure stability.  A very important point
is that the covariates in $\mathbf{X}(t)$ should be exogenous, which
excludes e.g. a temperature-driven heating power.  Note that with
$\alpha = 0$ and no covariates, the temperature difference $Z(t)$
would return to zero, which should make sense on physical grounds.
Now the corresponding approximate DT model for $Z_t$ at times $t=1$,
$2$, $\dots$ is
$$
  Z_t = \alpha + \beta Z_{t-1} + \boldsymbol{\gamma}^\top \mathbf{X}_t
  + \varepsilon_t,
$$
with $\beta := 1 - \theta$, so $\beta < 1$ should hold. 
This model slightly differs from your one because of the use of $Z_t$ in
place of $T_t$ at left hand side and, more importantly, because at the
right hand side we have the lagged $Z_{t-1}$ where you used
$Z_t$. This is an ARMAX model, with parameters $\alpha$, $\beta$, $\gamma_i$ and
$\sigma_{\varepsilon}^2$ which can be estimated by Maximum Likelihood
or Least Squares. In the second case, the parameters can be obtained
by a regression similar to the one you used.
To forecast $T_{t+1}$ using observations up to time $t$, you need
forecast $Z_{t+1}$ using a forecast of $\mathbf{X}_{t+1}$ 
$$
  \widehat{Z}_{t+1|t} = \alpha + \beta Z_t + \boldsymbol{\gamma}^\top \widehat{\mathbf{X}}_{t+1|t}
$$
where the parameters are replaced by their estimated value, and then
simply add the forecast of $Z_{t+1}$ to that of $T_{\textrm{a},t+1}$
according to $\widehat{T}_{t+1|t} = \widehat{T}_{\textrm{a},t+1|t} +
\widehat{Z}_{t+1|t}$. Multistep forecasts $\widehat{T}_{t+h|t}$ can be
derived as well for $h >1$.
Although general ARMAX functions exist in several R packages, they do
not all allow this simple parametric form. So it may be easier to keep
on writing your own estimation and forecasting functions using the
specific form of your model, unless you want to experiment with some
model changes. You can also try the R package named dynlm (on CRAN)
working on time series (ts) objects rather than on numeric vectors.
