Predictions for future n steps

[Navdeepuniyal]

I recently found this work and would like to use it for predicting the future 'n' values in a regression.
e.g. I have a dataset from t=0 to t-m and I want to predict the value of Y for steps t+n
Is it possible using NGBoost?

If yes, could someone please point me to an example?

[alejandroschuler]

Sure, this is possible, but it has more to do with how you shape your data than how to use ngboost. I'm assuming here that for a single unit of observation i in your dataset you have [yi0, yi1 .... yit ... yit+n] so you are trying to use the values [yi0 ... yit] to predict yit+n for some fixed n. At the moment there isn't a simple way to predict all of [yit+1 ... yit+n] in one shot, so you have to use one model for t+1, another model for t+2, and so on up to t+n (*footnote). So proceeding with some fixed n, what you need to do is shape your data so that you have a matrix X' and a vector Y'. Each row of X' will be X'i = [yi0, yi1 .... yit] and the corresponding element of Y' is Y'i = yit+n. Now simply feed X' and Y' into ngboost and you'll have a model capable of predicting yt+n from any vector [y0, y1 .... yt]. You can also include additional time-varying or static predictors as different rows of X' if you have them.

It's important to understand, however, that this model has absolutely no notion of "time". As far as it knows, y1, y2, etc. are unrelated numbers that may or may not have anything to do with each other and which may be combined in whatever way is best to predict yt+n. It also does not know or care that yt+n represents the same "measurement" as yt, just at a future time. Whether or not that is a feature or a downside of the approach depends on your goals and perspectives.

(*footnote) while there is not a currently-implemented approach that I would advise for simultaneous prediction of [yt+1 ... yt+n], it is absolutely possible in theory. As long as you can parametrize the vector [yt+1 ... yt+n] with finite-dimensional vector of parameters (e.g. with an AR1 model) then you can apply the ngboost algorithm to estimate them and therefore get probabilistic prediction over [yt+1 ... yt+n] simultaneously. Technically you could do this with the multivariate normal, but if n is large then that's a ton of parameters to fit and it will take forever. (and it ignores linear time). So probably best to think about an autoregressive distribution for this problem, which we don't have implemented and which I had actually not thought of until now.

[alejandroschuler]

(@MikeOMa might be interested in this conversation)

[MikeOMa]

I have thought about the ideas in the footnote. I was looking for low-rank positive semi-definite representations so you could do a larger output dimension with the multivariate normal without gaining that many parameters. My search didn't come up with anything particularly promising though.

The multivariate normal distribution would probably work well for n=2. It has n parameters for the mean which I think is unavoidable and (n^2+n) /2 for the covariance matrix. I haven't looked at any results for anything other than n=2 yet. I would be surprised if the multivariate normal did well for n>3 but who knows. There is an example of the multivariate normal in the examples directory of the github repository with n=2 if that is of interest. It is not a forecasting example but if you just make an array with vals[:,0] = y[t+1], vals[:,1] = y[t+2] the example would be similar.

I definitely agree that the benefit is in the multivariate prediction interval side! I'd say more often than not these correlation approaches will give worse metrics than the fitting n models strategy in 1 dimension. I am often wrong though.

Thoughts for the AR idea:
If what was wanted is a classical AR(1) covariance matrix on the outputs. That would be Cov(y_[t+k],y_[t+k]) = sigma^2 Correlation(y[t+k],y_[t+l]) = rho^abs(k-l). So just 2 covariance matrix parameters and n mean parameters.
I think there is a very nice equation for the Cholesky decomposition of the inverse of an AR(1) covariance matrix which should make the gradient and Fisher information calculations very straightforward. I had a quick look and I couldn't find it though. If anyone's interested in implementing it I'll do another search!

I think a AR type approach with unequal variance would be a bit more flexible. I was thinking something like Cov(y[t+k], y[t+k]) = sigma[k]^2, Corr(y[t+k],y_[t+l]) = rho^abs(k-l)
So you'd only have 2n+1 parameters to fit. I think the calculations for this would be slightly harder though. For n=2 this approach is completely equivalent to just using the current multivariate normal distribution implementation.

Edit:
The AR(1) idea was wrong! I was thinking of a loss function for the 1 step ahead errors. In the case of the k step ahead forecast you would probably want to structure the distribution based on the k step ahead AR forecast distribution. Even though there is an equal variance, the errors get larger with the forecast horizon. I am not sure if this would be straightforward.

[Navdeepuniyal]

Thank you for the insights, it is helpful.
I have another question related to the MultiOutputRegressor from sklearn. I am trying to predict multiple outputs using the same inputs. I see in MultiOutputRegressor, the predict_dist method is not supported. Is there any other way to achieve this?

[MikeOMa]

I would just make a wrapper class personally, edited an example from https://scikit-learn.org/stable/auto_examples/ensemble/plot_random_forest_regression_multioutput.html#sphx-glr-auto-examples-ensemble-plot-random-forest-regression-multioutput-py

# Author: Tim Head <betatim@gmail.com>
# Edited by: Michael O'Malley ,MikeOMa 
## Replaced the RandomForest multioutput with an ngboost example
# License: BSD 3 clause


import numpy as np
import matplotlib.pyplot as plt
import ngboost
from sklearn.base import clone
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.multioutput import MultiOutputRegressor
# Create a random dataset


rng = np.random.RandomState(1)
X = np.sort(200 * rng.rand(600, 1) - 100, axis=0)
y = np.array([np.pi * np.sin(X).ravel(), np.pi * np.cos(X).ravel()]).T
X_train, X_test, y_train, y_test = train_test_split(
    X, y, train_size=400, test_size=200, random_state=4)

max_depth = 30
regr_multirf = MultiOutputRegressor(RandomForestRegressor(n_estimators=100,
                                                          max_depth=max_depth,
                                                          random_state=0))
regr_multirf.fit(X_train, y_train)


class multioutput_ngboost:
    def __init__(self, estimator):
        self.estimator = estimator
        self.models = None
        self.n_models = None
    def fit(self, X,Y):
        """
        X: N x p array-like
        Y: N x d np.array
        """
        self.n_models = Y.shape[1]
        self.models = [clone(self.estimator).fit(X,Y[:,i]) for i in range(self.n_models)]

    
    def predict(self, X):
        """
        X: N x p arraylike
        Returns:
            a Nxd array containing predictions
        """
        return np.array([model.predict(X) for model in self.models]).T
    
    def pred_dist(self, X):
        """
        X: Nxp array-like
        Returns:
            List of ngboost distributions
        """
        return [model.pred_dist(X) for model in self.models]
    
    def score(self, X, Y):
        return np.mean([self.models[i].score(X,Y[:,i]) for i in range(self.n_models)])



estimator = ngboost.api.NGBRegressor(n_estimators=300, learning_rate=0.1)
ngboost_mo = multioutput_ngboost(estimator)
ngboost_mo.fit(X_train, y_train)

# Predict on new data
y_multirf = regr_multirf.predict(X_test)


y_ngboost = ngboost_mo.predict(X_test)
## Briefly check everything is working as it should
pred_dists = ngboost_mo.pred_dist(X_test)
assert all(pred_dists[1].mean() == y_ngboost[:,1])


plt.figure()
s = 50
a = 0.4
plt.scatter(y_test[:, 0], y_test[:, 1], edgecolor='k',
            c="navy", s=s, marker="s", alpha=a, label="Data")
plt.scatter(y_ngboost[:, 0], y_ngboost[:, 1], edgecolor='k',
            c="cornflowerblue", s=s, alpha=a,
            label="NGBoost multioutput %f" % ngboost_mo.score(X_test, y_test))
plt.scatter(y_multirf[:, 0], y_multirf[:, 1], edgecolor='k',
            c="c", s=s, marker="^", alpha=a,
            label="RF score=%.2f" % regr_multirf.score(X_test, y_test))
plt.xlim([-6, 6])
plt.ylim([-6, 6])
plt.xlabel("target 1")
plt.ylabel("target 2")
plt.title("Comparing a custom ngboost multioutput and random forest using sklearns multioutput")
plt.legend()
plt.show()

Is this roughly what you were looking for?

[Navdeepuniyal]

This is great. Thanks a ton @MikeOMa