Estimate local slope of a fitted model's response surface with respect to time

Fits a user-supplied model to a univariate time series (using time index and an auxiliary random covariate as predictors), then numerically differentiates the fitted model's predictions with respect to time using a central finite-difference (perturbation) scheme. Slopes are computed for every historical time point as well as for h future time points, making this useful for characterizing local trend/momentum in the series both in-sample and over a forecast horizon.

estimate_theta_slope(fit_func, predict_func, y, h = 0, ...)

Arguments

fit_func

A function used to fit the model. It is tried successively with a formula interface (y ~ ., data = ...), an x/y interface, and a positional (x, y) interface, so it should support at least one of these calling conventions (e.g. lm, randomForest, glmnet, etc.).

predict_func

A function used to generate predictions from the object returned by fit_func. It is called as predict_func(object, newdata) and, if that fails, as predict_func(object, as.matrix(newdata)). Must return a numeric vector of predictions.

y

A numeric vector containing the univariate time series of historical observations.

h

Integer; forecasting horizon, i.e. the number of future time points (beyond the length of y) for which slopes should also be computed. Defaults to 0 (in-sample slopes only).

...

Additional arguments passed on to fit_func.

Value

A numeric vector of length length(y) + h containing the estimated local slope (numerical derivative of the fitted model's predictions with respect to time) at each historical and future time point, in chronological order.

Details

The function builds a data frame with columns y (the response, padded with NA for the h forecast periods), t (an integer time index running from 0 to n + h - 1), and z (an i.i.d. standard normal covariate used only to satisfy models that require more than one predictor). The model is fit on the first n rows (the historical data) via fit_func, trying in order:

  • the formula interface fit_func(y ~ ., data = df_train, ...);

  • the matrix/vector interface fit_func(x = cbind(t, z), y = y, ...);

  • the positional interface fit_func(cbind(t, z), y, ...).

If none of these succeed, the function stops with an error.

For every time point (historical and future), an adaptive step size h_eps = max(zero^(1/3) * |t|, zero) (with zero = 1e-4) is used to perturb t up and down while holding z fixed. The fitted model is used to predict the response at t + h_eps and t - h_eps, and the central-difference slope (fx_plus - fx_minus) / (2 * h_eps) is returned for each time point.

A fixed random seed (123) is set internally so that the random covariate z — and therefore the resulting slopes — are reproducible across calls.

Examples

if (FALSE) { # \dontrun{
set.seed(1)
y <- cumsum(rnorm(50))
slopes <- estimate_theta_slope(
  fit_func = lm,
  predict_func = predict,
  y = y,
  h = 10
)
plot(slopes, type = "l")
} # }