Nixtla's Open Source Time Series Ecosystem.

Nixtla Community

That's the expected behavior of `falback_model.`
 :) 
Please let us know if something fails. :raised_hands:

image.png

Hello everybody. I am testing Nixtla for some hierarchical forecasts and I am getting the following error. Any one know what I can do? Thanks

hey <@U04836R2NGZ>! I’m suspecting that the problem is related to `MinTrace(method='mint_shrink')` . Could you please try removing it to know if that’s the method that breaks the pipeline? :raised_hands:

Hello Fede. Worked. Thank you. Also, I just recently started using Nixtla. What a great job you guys have done. :+1::skin-tone-4: congrats.

Thank you, <@U04836R2NGZ>!  I have a good idea of what might be happening with the `mint_shrink` bug. We will fix it soon.

Please let us know if you have any other issues or special features you would like `hirarchicalforecast`  to have. :slightly_smiling_face:

Hi <@U04836R2NGZ>,

I am working on debugging the MinTrace-shrinkage method.
I wonder if the data that you are using for your experiments is public, that would help me to explore the numerical instability causes.

Hello Kin, I can send it to you. No problem. But I can send it tomorrow because I’ve left the office (8 hour time diff.)

That sounds good, blind debugging is on steroids :+1: 

Hi <@U03UANPFB8U> and <@U04836R2NGZ>

I added numerical protections to `MinTrace` , available now in <https://github.com/Nixtla/hierarchicalforecast|HierarchicalForecast's main>.

The origin of the error are series with very low insample residual variance. Example `y_model=[1,1,1,1] y_hat=[1,1,1,1]`. You can force MinTrace to run in such scenarios with the following:

• I added the `mint_shr_ridge` parameter to `MinTrace` class, that forces the covariance matrix to be invertible. In case you want to numerically force the method.
• Another trick that helped me in the past on these ill conditioned variance matrix is to add white noise to `Y_df['y']` that works almost like the `mint_shr_ridge` , but in a stochastic version.
Hope this helps

Thank you <@U038ZT9UDGA> for the fix and update!

Thanks for reporting the error. Let us have a look and come back to you. 

Hi <@U03UANPFB8U>; the `HierarchicalLoad` function is meant for benchmarking or research purposes. The idea is to have a function that makes it easier to load famous datasets like _tourism_ or _Labour_.

If you already have forecasts at all different levels, you just need the `S` matrix and the `tags` to use any reconciliation method.
The `aggregate` function assumes you just have the lowest-level hierarchy. Currently, we don’t have an explicit function to generate the `S` Matrix or `tags` for the case where the user already has forecasts for every level.

However, you can do the following “hack” and it should work fine.

1. Create a `Y_fitted_df_bottom` with just the bottom-level ids. I.e, filter your current `Y_fitted_df` so it has just the ids from the lowest level.
2. Use the aggregate function to generate your `S`  Matrix and `tags`. (You can trash the resulting df.) Something like:  `Y_fitted_df_trash , S_fitted, tags_fitted = aggregate(Y_fitted_df_bottom, hiers)`
3. Use your original `Y_fitted_df,` and the newly generated `S` and `tags` to call the <https://nixtla.github.io/hierarchicalforecast/methods.html|reconciliation method >of your choosing.


Hi <@U03UANPFB8U>

This is an example where I did the hack to recover the `S_df`
from a complete hierarchy, you would need to adapt this example to your needs:
```# Reading Data
Y_hat_df = pd.read_csv('Y_hat_df.csv')
Y_fitted_df = pd.read_csv('Y_fitted_df.csv')

# Creating bottom df from Y_fitted
unique_ids = Y_hat_df.unique_id.unique()
level_count = np.array([u_id.count('/') for u_id in unique_ids])
is_bottom = (level_count == max(level_count))
bottom_level =  np.array([u_id[4:] for u_id in unique_ids])
bottom_df = pd.DataFrame(dict(unique_id=unique_ids,
                              is_bottom=is_bottom,
                              total_level=['SAU']*len(unique_ids),
                              bottom_level=bottom_level))
bottom_df = bottom_df[bottom_df.is_bottom==True]
bottom_df['ds'] = 1
bottom_df['y'] = 1
bottom_df = bottom_df[['unique_id', 'total_level', 'bottom_level', 'ds', 'y']]

hierarchy_levels = [['total_level'], ['total_level', 'bottom_level']]
_, S_df, tags = aggregate(df=bottom_df, spec=hierarchy_levels)```
This is a nice idea for a functionality of the library.

Hi <@U04816WA7N3>

If you are not using the NeuralForecast wrapper and only the PyTorch model these lines do the trick:

```from pytorch_lightning.callbacks import TQDMProgressBar
from neuralforecast.tsdataset import TimeSeriesDataset, TimeSeriesLoader

temporal_df = data['temporal'].copy()
temporal_df = temporal_df.groupby('unique_id').head(228).reset_index()
dataset, *_ = TimeSeriesDataset.from_df(df=temporal_df, sort_df=False)

model = GMM_TFT(input_size=12*4,
                hidden_size=64,
                h=12,
                K=2,
                k_cont_cols=['distance_month', 'y_[lag12]'],
                k_cont_inp_size=2,
                batch_size=32, #4,
                windows_batch_size=256,
                check_val_every_n_epoch=1,
                callbacks=TQDMProgressBar(refresh_rate=1, process_position=0),
                loss=GMM(quantiles=np.arange(0,100,5)[1:]/100),
                #loss=GMM(quantiles=np.arange(0,100,2)[1:]/100),
                #learning_rate=0.1,
                learning_rate=5e-3,
                max_epochs=2,
                enable_progress_bar=True)

# Fit and predict
model.fit(dataset, val_size=12, test_size=12)
Y_hat = model.predict(dataset=dataset)
Y_hat = Y_hat.reshape(555, 12, len(model.loss.quantiles))```

We just added static variables to the model and homogenized the inputs, here is:
• <https://nixtla.github.io/neuralforecast/models.gmm_tft.html|GMMTFT documentation>.
• <https://github.com/Nixtla/neuralforecast/blob/main/neuralforecast/models/gmm_tft.py|GMMTFT code link>.
Let me know how it goes.

Here is a link to <https://pytorch-lightning.readthedocs.io/en/stable/api/pytorch_lightning.callbacks.EarlyStopping.html|early stopping documentation>.
you would need to send the callbacks with early stopping like a list.

Keep in mind that good performance in validation does not correlate as well in M5. If you are working with that dataset.

Thank you! I have used pytorch-lightning early stopping callback before - if the syntax is standard, I shouldn't have any issues. As far as my dataset - I'm using a simpler dataset with a continuous valued target (drug prices), since that seems more appropriate for the GMM loss. I think validation performance ought to be relevant but we'll see!

<@U04816WA7N3>, I would recommend you to first try <https://nixtla.github.io/neuralforecast/examples/getting_started_with_nbeats_and_nhits.html|NHITS/NBEATS with MAE>,
then move towards MQLoss and then to GMMTFT, that is still research code.

Thanks, I'll do that! Also - it seems the GMM_TFT class needs to be updated with the new parameters for TFT (tested with the documentation AirPassengers example):

Thanks <@U04816WA7N3>
I created this <https://github.com/Nixtla/neuralforecast/issues/308|github issue>, will get to it in the week

Is this the full code? I can’t see what would cause this immediately

not really, what I am trying to do is, since I have a large of models to run and apply, I want to break model training and forecasting into 2 processes and use 2 Pandas UDF to make it happen. But I am not sure how can I log the statsforecast model in ML Flow and how can I call them.

I’ll make an example for you later today

here is the code i have right now just in case. the structure was built based on an reference I had from sklearn, but it seems that log_model function is not applicable here.

Ok so I have an end to end script that is working: <https://github.com/kvnkho/demos/blob/main/fugue/fugue-nixtla-mlflow.py>

I don’t know MLFlow well enough to get the nested runs working yet. I can probably revisit it more Wednesday, but I think this should give a clear picture of where this is going. Happy to hop on a call and discuss. I’m just a bit busy today and tom. Also, I am not with the Nixtla team if ever they correct me lol, am just a friend. I am part of the <https://github.com/fugue-project/fugue/|Fugue> team where we work on bringing Python and Pandas code to Spark easily (as you’ll see in the end of the script). If you don’t want to use Fugue, it’s fine, you can still use something like pandas_udf. It will just be some extra code to re-write (especially around Schema).

I am presenting in PyData NYC next week and this is a good use case so I’ll be use to understand it more and get it working by then. If you do get the nested client working, I’d love to see it!

Thank you so much for sharing Kevin, I will look into the code and fit the same on my dataset to test. Appreciate you efforts with the script!!

<@U03SA41952R> Hi Kevin, I spent some time to look into the code and have some follow-up questions if you don't mind. I am trying to break training models and predicting as 2 process, since when they are combined as one, the computation time increase drastically. I learned this from a different project. With this situation, (1) how can I use the saved model.pkl file and apply it for forecasting? I need  exogeous variable in this model, thus I am planning to have model.fit() in the 1st process, and model.predict(X_future) in the 2nd process. (2) If pickle file cannot be used in this way, is there any other method that can help me achive my goal here? I know that sklearn has load_model(uri) but not sure if there is anything similar in statsforecast?

Hey! Responding to general chat so Nixtla team can chime in. So my understanding with `statsforecast` is you can either use `.forecast()` which is a straight forecasing or `.fit()` and `.predict()` . If you want to decouple training and prediction time, just call the `.fit()` method and then you can save that pickle object to be used later.

The pickle should be able to be used with exogenous regressions because it saves the model weights and info. That is how it was done with <https://scikit-learn.org/stable/model_persistence.html#python-specific-serialization|sklearn before also>. I don’t believe `statsforecast` has it built in yet but you can just make your own saving and loading functions pretty easily.

The <https://nixtla.github.io/statsforecast/examples/exogenous.html#sklearn-syntax|exogenous docs> say that the class can handle the future exogenous regressors.

On local, the `mlflow.sklearn.load_model` just uses <https://github.com/mlflow/mlflow/blob/master/mlflow/sklearn/__init__.py#L425-L453|pickle> under the hood as well :slightly_smiling_face:

this might be another dumb question...but how can I get the location of the pickel file logged in the experiment? like, when I do pickle.load() what would be the location of model_filename?

I don’t know MLFlow enough but I think you get the MLFlow client, supply experiment name, get the artifact location, make the URL and then load the model

Ehh it doesn’t look too good. Seems like you need to know a bit about where the path exists: <https://stackoverflow.com/a/72625756/11163214>

Also see the method signature <https://www.mlflow.org/docs/latest/python_api/mlflow.artifacts.html#mlflow.artifacts.download_artifacts|here>

yep I saw the same post and got more confused :smile: will try this load_dict method! thanks a ton!

No not the `load_dict` . The one above `download_artifacts`.

I am so confused whether to use the Client or not. Last time I used MLFlow heavily was 3 years ago but we didn’t reply on it to manage the paths. We just using it to track experiments but for keeping track of paths of saved models we always had a “latest_model” path that was constant and we just always pulled from there.

hey <@U0423GYCJSY>! With `StatsForecast`  you can recover the fitted values using `fitted=True`  in the `forecast`  method. After that, you just have to call `forecast_fitted_values()`  without parameters to obtain the insample forecasts. Here’s a complete example,

```from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA, ETS, RandomWalkWithDrift
from statsforecast.utils import AirPassengersDF as Y_df

sf = StatsForecast(df=Y_df, 
                   models=[AutoARIMA(season_length=12), ETS(season_length=12), RandomWalkWithDrift()],
                   freq='M')

forecasts = sf.forecast(h=12, fitted=True)
insample = sf.forecast_fitted_values()```

Thanks for the quick response Fede! I will test the same code on my end :slightly_smiling_face:

just a side question, for ETS, we cannot get confidence intervals like level=90, is that correct?

yes, that’s correct, currently we don’t have levels implemented for ETS but we are working on including them

hmm I tested the code you shared above and it's still giving me the error like this..

ahh ok, are you also passing `level`  to `forecast`  right?

yeah...is this the root cause of the error?

yes, since some models dont have levels implemented there is an error trying to construct the dataframe of insample forecasts. But, ideally statsforecast should return the levels only for models that have them implemented

hmmm I removed the level argument but the error persists...

oh that’s unexpected, may I see your code?

wow I run the code again and it worked thru this time!! thanks for the help!

nice! i will fix the levels error and i’ll let you know once it is resolved

that would be awesome! please keep me posted :slightly_smiling_face:

<@U0316L4HWQ7> Hi Fede, a follow-up question on the thread above, if one model among the 3 (ARIMA, RWS, ETS) cannot get a model from the training data, it's it possible to pass some random value in the result dataset and continue running, instead of shutting the run down and throwing error?

I am having this question, since I want to find the best model for a lot of datasets, if one model is throwing an error, this certain dataset won't have any other models available.

hey <@U0423GYCJSY>! StatsForecast includes an argument called `fallback_model.`

The main idea behind it is that you can pass a safe model in case that a main model fails. So, StatsForecast tries to fit the main models for each time series and if something fails, it will fit a safer model, for example, `Naive().`

The code would be,

```StatsForecast(models=[AutoARIMA(), ETS()], fallback_model=Naive(), extra params) ```
Please let us know if that works for your use case. :slightly_smiling_face:

This is an awesome function!! Thank you for sharing Fede!

Is it averaged for all timeseries or just one?

Hey <@U03SA41952R>!

The `fit` method stores a `fitted_` attribute containing the fitted models for each time series. Each fitted model is different. For example, the naive model stores the last time series observation, while Arima stores other objects required for forecasting. Currently, we do not have the functionality to use a different time series for forecasting. The `predict` method only receives the horizon to forecast, so if the fitted model is saved, the forecasts produced by it will be the same each time `predict` is called. We are working on adding the functionality to use new data. For the moment, we recommend using the forecast method.

Ah how does it work then for exogenous variables? What if `fit()` has one? How is it supplied to predict?

Ah I see, you just supply it also. Ok got it! Thanks!

combined.parquet

```import pandas as pd
from statsforecast.distributed.utils import forecast
from statsforecast.distributed.fugue import FugueBackend
from statsforecast.models import AutoARIMA

df = pd.read_parquet(f"{WORKING_DIR}/combined.parquet")[["unique_id", "ds", "y","sell_price"]]
result = forecast(df, 
                  models=[AutoARIMA(season_length=7)], 
                  freq="D", 
                  h=7,
                  parallel=FugueBackend("dask"))
result.compute()```

Maybe I have something wrong but this is pretty close to the tutorial and I get the same error
```from statsforecast import StatsForecast

df = pd.read_parquet(f"{WORKING_DIR}/combined.parquet")
df = df[["unique_id", "ds", "y","sell_price"]]

model = StatsForecast(
    df=df,
    models=[AutoARIMA(season_length=7)], 
    freq='D', 
    n_jobs=-1
)
model.forecast(7)```

hey <@U03SA41952R>! Future variables must be included in the `forecast` method to use exogenous variables. That is, the code would be:

```model.forecast(7, X_df=future_exogenous) ```
For example, if you are using weather as an exogenous variable, you must include the weather for the days you generate the forecasts.

Here’s an example in a local setting: <https://nixtla.github.io/statsforecast/examples/exogenous.html>

For a distributed setting, we’ve opened a PR to add the functionality to fugue’s backend: <https://github.com/Nixtla/statsforecast/pull/300>

exogenous.html

Hello, I'm getting this error for a very simple dataset

I checked that there are no nulls in my data set

I think the problem is that you have a unique_id for every timestamp. Try changing the unique_id value to something like “1” or “Series 1" if you are just using 1 series.

StatsForecast can also handle multiple unique_ids automatically. If that is what you are looking for.

Also, if your data is weekly and you have more than 1 year of obserbations you can set the `season_length`  parameter of the AutoArima to 52.

So, it would look something like this

```df_agg2 = df_agg2.assign(unique_id='1')

sf = StatsForecast(df=df_agg2,
                   models=[AutoARIMA(season_lenght=52)],
                   freq='W',
                  n_jobs = -1)
forecasts = sf.forecast(7)
forecasts.head()```

Ah yeah the unique_id all being set to 1 worked

What if the seasonality is more monthly than weekly though?

For setting the optimal season_lenght you can check this article: <https://robjhyndman.com/hyndsight/seasonal-periods/>

You are welcome. Also, feel free to try more models and see what works best.
You just have to change this line:
```models=[AutoARIMA(season_lenght=52), ETS(season_lenght=52), ... ],```


Happy to help if you need anything else :slightly_smiling_face:

If I understand the question correctly, the second option is more efficient. 

I'll give it a try and come back if something seems off. Thanks for the quick reply!

Hello again. So I have a bit of an issue with using both the .forecast and .cross_validation methods in testing. In both cases, it seems to be fitting the model to the new data, which during testing might not be desirable, nor efficient. Is there a way to force it to use the trained parameters instead of retraining? Am I missing something in my approach?

ChainedWindows.gif

(Credit to *joaquinAmatRodrigo* for the second image)

Currently, the cross-validation class does the first thing.

If I understand correctly, you want to train a model on a subset of the data and forecast the future without retraining on newly available data. Maybe to save some computing resources and maybe because the model’s accuracy will not decrease so much. Right?

If that’s the case, that functionality is currently not available. But we are working on it.

However, in our experience, it makes sense to retrain “every time” you have new relevant data. Since statsforecast is highly efficient, even thousands of series take a couple of minutes. If you are using some of the Auto models and want to improve the training speed further, you can restrict the number of models. For example, `ETS (season_length:int=1, model:str='ANN')` will explore only a simple exponential smoothing.

Please also check this issue and like it if is are relevant to your use case: <https://github.com/Nixtla/statsforecast/issues/287>.

That's exactly it. So cv is doing like the first image, but with overlapping windows (unless one specifies step_size=horizon, then we get exactly the first image), right? With the overlapping windows then, if I had a test et of 100 points, it would effectively train 100 new models. This test would give me an estimate of what my model will be like in production if I retrain it every single day. Is that accurate?

I'm glad to hear you are working on the backtesting without refit. It would allow one to also measure how often the models need to be retrained or if the trained model is relatively stable through time. I'll be keeping an eye out for that functionality.

Thank you also for the tip on making the training more efficient. Is there anything else I should keep in mind for speed-up if I start doing a big number of refits. currently I only use num_jobs=-1. Should I think of installing any external libraries, do I need to setup use of GPU with neuralforecats models or is this all out of the box already?

&gt; That’s exactly it. So cv is doing like the first image, but with overlapping windows (unless one specifies step_size=horizon, then we get exactly the first image), right? With the overlapping windows then, if I had a test et of 100 points, it would effectively train 100 new models. This test would give me an estimate of what my model will be like in production if I retrain it every single day. Is that accurate? (edited)
Yes!

We did a new release in `HierarchicalForecast` , show some love community :heart:

<https://twitter.com/fede_gr/status/1592621393975586817>

TCN work looks great btw, tossing this here in case there’s interest in it: <https://lnkd.in/eBcGjh3V|https://lnkd.in/eBcGjh3V>

We played around with neural ODEs a bit earlier this year but found them very slow to train and quite unstable, really excited about the new research direction!

Looks like it’s implemented in the `AutoARIMA` class, but still getting fairly low predictions: `AutoARIMA(stepwise=True,season_length=13,blambda=0.4,biasadj=0.4)`

Hi <@U03UTAG6GNP>, thanks for the question. We will have a look and come back to you with an answer.

The `predict` and `forecast` method takes two arguments: forecasts the next `h` (for horizon) and`level`.


* `level` (list of floats): this optional parameter is used for probabilistic forecasting. Set the `level` (or confidence percentile) of your prediction interval. For example, `level=[90]` means that the model expects the real value to be inside that interval 90% of the times.

The forecast object here is a new data frame that includes a column with the name of the model and the y hat values, as well as columns for the uncertainty intervals.

Please check this guide on <https://nixtla.github.io/statsforecast/examples/uncertaintyintervals.html|Arima and Predictions Intervals> and let us know if something doest not work.

fcsts_df = fcst.forecast(h=12, fitted=True, level=(90, 10))

Here's a tutorial for StatsForecast's AutoARIMA <https://nixtla.github.io/statsforecast/examples/getting_started_with_auto_arima_and_ets.html>

Now, if you want something more theoretical, Hydman's book explains how AutoARIMA works. <https://otexts.com/fpp2/arima.html>

as197 1.png

Hi <@U0471GWET5F>, thanks for posting here. We really enjoyed your talk at the dev days of sktime and felt humbled that you mentioned us :) 

Just yesterday we were discussing about the mentioned recursion method and we’re intrigued by it. 

As we understand it, the massive speed gains apply only to certain ARIMAs and came to the conclusion that at the moment we are not sure about the trade off between implementation effort vs those gains. That being said, we are very happy to discuss the issue with you. You opinión is going to be surely very helpful. (Have you played around with it or have some intuitions?)

Regarding bigger season lengths we have experienced good results using MSTL (with ARIMA). 

Thanks for the response Max! I've done a small benchmark (~50 time series) comparing Gretl (the C econometrics library)'s auto_arima with AS 197 against statsforecast and found the following: the variance in fit times is much lower (between 1-10 seconds vs 3-60). Once again, it seems that statsforecast's auto_arima if "fast enough" if we know beforehand that the time series is relatively well-behaved (e.g. AR(1)MA(1)) or (as you suggested) applying MSTL decomposition beforehand. That being said, no other implementation in Python compares!

To clarify, I found that the longer fit times come from more "complex" seasonal time series, i.e. when the step-wise search has to look at higher orders before early-stopping. So in production, I've lowered the max p, q, P, Q by half. <https://www.youtube.com/watch?v=uyIlAO390v4|https://www.youtube.com/watch?v=uyIlAO390v4>

Also <@U035DR8HD6D> thanks for watching the talk and looking into the suggestion! Hopefully I've made a convincing case to move towards global models in a vertically scalable setting. Looking forward to seeing how `neuralforecast` develops!

P.S. Your MSTL-ARIMA suggestion(ala <https://www.census.gov/data/software/x13as.html|https://www.census.gov/data/software/x13as.html>) should probably be the standard benchmarking "AR" technique over standard auto_arima....it is just 2-3 factors slower than all other benchmarking models.

Do you think a thin wrapper around such a model such as x13as is a good idea? 

P.S.S. Have you seen this? <https://www.prnewswire.com/news-releases/arima-on-steroids-anodot-gets-new-us-patent-for-lightning-fast-anomaly-detection-solution-301589990.html>

I really enjoyed reading this thread!
About Anodot, I’m curious about what they are doing.

&gt; Do you think a thin wrapper around such a x13as model is a good idea?
Great idea :slightly_smiling_face:

Hi! For production settings we recommend using the forecast method instead of the fit predict ones. 

Also, you should explore the possibility of using other models like MSTL or ETS. 
This newly released guide might be helpful. 
<https://nixtla.github.io/statsforecast/examples/getting_started_complete.html|https://nixtla.github.io/statsforecast/examples/getting_started_complete.html>

<@U035DR8HD6D> <https://github.com/Nixtla/statsforecast/tree/main/experiments/arima_prophet_adapter> -&gt; I see you here using fit and then make future dataframe and then predict. Could you be more specific what you mean by using forecast?

Secondly, we want to test, on our dataset, how much of a speed and accuracy improvement we get

Hi Andrei, I will try to give you a complete answer soon.

Do you have any code that you could share so i can be more specific?

Yes, I will ask the guys from our team to  send over the code. But most plainly put. compared directly Autoarimaprophet with prophet, with uncert intervals = 0 and some fut regressors. The dataset is weekly and kinda small.

Is there a specific reason you want to compare Prophet and AutoArima? I would highly suggest finding appropriate models for each combination of series (e.g. product and store).

For some series, sparse models like ADIDA or Croston might be very successful.

we have weekly data and it`s not sparse. And we want to use it because 1. we have seasonality and trend + maybe some holidy effect (product demand forecasting). we have some future regressors and  dataset is not that large to use deep learning models

For that use case, I suggest comparing Prophet against MSTL.
MSTL decomposes seasonality and trends but can fit multiple seasonalities. It’s also faster than ARIMA.

That’s good for your use case because with weekly data, you might find monthly patterns (last week of the month) as well as yearly patterns (last week of December)

Would you mind looking at <https://nixtla.github.io/statsforecast/examples/multipleseasonalities.html|this guide> and telling me if that’s something you would try?

Hi everyone, I joined this Nixtla community slack channel just now. I am working on hierarchical time series project. I am following this <https://www.kaggle.com/code/konradb/ts-8-hierarchical-time-series/notebook> to implement the hierarchical model.

I am using the following command to install the required libraries.
*pip3 install hierarchicalforecast statsforecast datasetsforecast*

Hi <@U04CVB9M7L1>! I’m sorry to hear and happy to help.  Could you share the error messages and maybe some code so we can help you out ?

I tried to install hierarchicalforecast library using the following command.
*<https://github.com/carlomazzaferro/scikit-hts-examples/blob/master/notebooks/visnights.ipynb|pip3 install hierarchicalforecast>*